Atelier R ANF FOSS : La régression géographiquement Pondérée ou GWR

Thierry Feuillet (UMR IDEES, Université de Caen-Normandie)

Grégoire Le Campion (UMR Passages, CNRS)

Objectif de l’atelier : Être en mesure de réaliser et interpréter une régression géographiquement pondérée (GWR).

Prérequis : Connaissances de base en analyse de données statistiques, avec au moins une compréhension de ce qu’est une corrélation et une régression.

1 Pourquoi la GWR ?

Lorsque l’on souhaite dépasser la simple caractérisation d’attributs liés à des individus statistiques, on fait appel à des méthodes de modélisation explicitant des relations statistiques. La méthode la plus employée, principalement en sciences sociales, pour mesurer et analyser la nature des relations et des effets entre deux ou plusieurs variables, est le modèle de régression linéaire.

Le principe de la régression linéaire est de modéliser la variable que nous souhaitons étudier (aussi appeler variable dépendante, VD) comme une fonction linéaire des variables que nous aurons définies comme explicatives de la VD (aussi appelées variables indépendantes, VI). Lorsque l’on s’intéresse à un phénomène social avec une emprise sur un espace, la régression linéaire pose plusieurs problèmes :

Le premier est empirique. La régression linéaire nous permet d’obtenir des coefficients (appelés betas β) et des résidus (notés epsilon ε). Ces β représentent l’effet de nos VI sur notre VD. Ces β sont considérés comme globaux, sans variation. Autrement dit, les modèles de régression linéaires considèrent que les VI interviennent de la même manière et avec la même importance sur l’ensemble de notre jeu de données. Si cette hypothèse peut être validée sur des populations statistiques définies aléatoirement et sans effet de structure a priori des VI ou de la VD, elle n’est que rarement vérifiée sur des données spatiales. En effet, les caractéristiques propres de chaque territoire (l’unicité de chaque lieu) impliquent que l’effet constaté en un lieu n’est pas forcément valable en un autre lieu de l’espace.

Concernant les prix des valeurs foncières que nous détaillerons par la suite dans cette fiche, on peut comprendre que la proximité au littoral, très prégnante en certains points de l’espace, ne joue absolument aucun rôle dans d’autres lieux. De même, certaines caractérisations du monde rural n’interviennent plus lorsqu’on se situe dans des milieux fortement urbanisés. Ainsi, les données spatialisées sont soumises à l’hétérogénéité spatiale : l’effet de nos VI va varier en fonction de l’espace. Un coefficient qui serait global et uniforme pour mesurer un effet paraît plus simple et donc tentant, mais non pertinent en géographie ; sur ce point nous pouvons nous référer à l’article de Brunsdon, Fotheringham et Charlton (Brunsdon et al. 1996). Ce concept d’hétérogénéité dans l’espace se traduit en statistique par celui de non stationnarité.

Le deuxième problème est statistique : chaque méthode statistique doit répondre à un certain nombre de conditions de validité. La régression linéaire ne fait pas exception. Trois conditions doivent être validées pour qu’une régression linéaire puisse être effectuée sans que l’interprétation des résultats conduisent à des raisonnements fallacieux :

Les individus statistiques doivent être indépendants
Les résidus doivent suivre une distribution normale
Il ne peut pas y avoir plus de VI que d’individus statistiques

Si les deux dernières conditions ne trouvent pas de matérialisation spécifique sur des données spatiales, la première quant à elle concrétise un problème récurrent sur les données en géographie. Par leur nature même, les données spatiales ne peuvent pas remplir cette condition fondamentale pour une régression classique. La première loi de la géographie de Tobler : “everything is related to everything else, but near things are more related than distant things” en est une traduction tout à fait parlante.

Le troisième problème est lié aux effets de contexte : on ne peut pas étudier des données spatiales sans considérer que les individus statistiques (les objets spatiaux) appartiennent eux-mêmes à des agrégats qui ont une influence sur la variable à expliquer. Ainsi, certains attributs de l’agrégat vont avoir une influence sur l’entité spatiale de cet agrégat.

Le quatrième problème est lié à la problématique du MAUP (Modifiable Area Unit Problem). Il concerne un problème d’échelle d’application de la régression et peut conduire à des observations erronées. Une corrélation constatée à une échelle peut être uniquement liée à l’agrégation réalisée à cette échelle, mais s’avérer erronée à une échelle plus fine, invalidant de fait la relation entre les phénomènes étudiés (Mathian, Piron 2001). Certaines agrégations peuvent également varier et cacher des relations entre individus (Bailey, Gatrell 1995).

La GWR ne répond pas à l’ensemble de ces problèmes mais va nous permettre de résoudre les deux premiers en intégrant la dimension spatiale de nos données tout en tenant compte de l’hétérogénéité (ou non stationnarité) de leur effet.

1.1 Les packages

Voici les packages que nous utiliserons :

# Chargement, visualisation et manipulation de la données
library(here)
library(DT)
library(dplyr)

# Analyse et représentation statistique
library(car)
library(correlation)
library(corrplot)
library(ggplot2)
library(gtsummary)
library(GGally)
library(plotly)

# Manipulation et représentation de la données spatiales (cartographie)
library(sf)
library(mapsf)
library(rgeoda) #permet en plus de calculer les indices d'auto-corrélation spatiale
library(RColorBrewer)

# Calcul du voisinage et réalisation de la GWR
library(spdep)
library(GWmodel)

1.2 Présentation et préparation des données

Nous avons cherché à traiter ici une variable présentant des caractéristiques spatiales fortes et qui rencontrent les deux problèmes exposés précédemment d’indépendance statistique et de non-stationnarité. Les prix de l’immobilier en France sont effectivement soumis à ces problèmes et peuvent être expliqués, au moins partiellement, par des variables quantitatives issues de données de l’INSEE. Nous avons traité l’information liée au prix de l’immobilier à l’échelle de Etablissements Publics de Coopération Intercommunale (EPCI) pour nous assurer un nombre d’observation suffisant dans chaque entité spatiale.

Les données du prix de l’immobilier par EPCI (prix médian au m²) sont issues des ventes observées sur l’année 2018, extraites depuis la base de données des notaires de France par Frédéric Audard et Alice Ferrari. Ce fichier a été simplifié pour ne conserver que les variables d’intérêts parmi une cinquantaine
Les données statistiques proviennent de l’INSEE (année 2019) : 9 variables ont été choisies pour leur potentialité à expliquer les variations des prix de l’immobilier, concernant la population, le logement et les revenus et niveaux de vie. Elles sont détaillées un peu plus bas.
Comme indiqué en introduction les données spatiales proviennent de la base ADMIN-EXPRESS de l’IGN en accès libre. La couche est utilisée est celle des EPCI de la base ADMIN-EXPRESS-COG édition 2022 par territoire pour la France métropolitaine.

Les données statistiques et du prix de l’immobilier ont été regroupées dans un même fichier CSV donnees_standr.csv.

Les données spatiales (géométrie des EPCI) sont au format shapefile dans le fichier EPCI.shp.

1.2.1 Chargement des données sur le prix de l’immobilier par EPCI

# On situe le dossier dans lequel se trouve nos données
csv_path <- here("data", "donnees_standr.csv")
# lecture du CSV dans un dataframe
immo_df <- read.csv2(csv_path)
# Pour visualiser les 10 1ères lignes
datatable(head(immo_df, 10))

Ce fichier est composé des 10 variables suivantes :

SIREN : code SIREN de l’EPCI
prix_med : prix médian par EPCI à la vente (au m²)
perc_log_vac : % logements vacants
perc_maison : % maisons
perc_tiny_log : % petits logements (surface < ?)
dens_pop : densité de population (nb habitants / km² ?)
med_niveau_vis : médiane du niveau de vie
part_log_suroccup : % logements suroccupés
part_agri_nb_emploi : % agriculteurs
part_cadre_profintellec_nbemploi : % cadres et professions intellectuelles

La variable SIREN nous servira de “clé” pour joindre ces données statistiques aux données spatiales, la variable prix_med sera la variable que nous chercherons à expliquer (VD), et toutes les autres seront nos variables explicatives (VI).

1.2.2 Chargement des données géographiques : les EPCI de France métropolitaine

Ces données proviennent de la base ADMIN-EPXRESS-COG de l’IGN, édition 2022. Le format d’entrée est le shapefile mais nous passerons par une conversion au format sf, ce qui nous permet d’utiliser le package mapsf, pour les prévisualiser :

# lecture du shapefile en entrée dans un objet sf
shp_path <- here("data", "EPCI.shp")
epci_sf <- st_read(shp_path)
## Reading layer `EPCI' from data source 
##   `/media/glecampion/Data/gregoire/atelierR_foss_gwr/data/EPCI.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1242 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 124110 ymin: 6049642 xmax: 1242428 ymax: 7110453
## Projected CRS: RGF93 v1 / Lambert-93
# visualisation des données géographiques
mf_map(x = epci_sf)

# et la table attributaire correspondante
datatable(head(epci_sf, 5))

1.2.3 Jointure des données géographiques et tabulaires

Les lignes vides ne nous intéressent pas, nous ne les conserverons pas car elles pourraient poser problème lors de la réalisation de nos analyses. On fait la jointure en ne gardant que les EPCI ayant une correspondance dans le tableau de données :

data_immo <- merge(x = epci_sf, y = immo_df, by.x = "CODE_SIREN", by.y = "SIREN")
nrow(data_immo)

## [1] 1223

Pourquoi certains EPCI n’ont-ils pas de correspondance dans notre tableau de données ? Il s’agit notamment des EPCI de la petite couronne de Paris, qui se superposent à l’EPCI de la métropole du grand Paris (les données sont donc bien présentes pour cette zone). Pour le reste, il s’agit de nouveaux EPCI ou bien d’EPCI ayant évolué et donc changé d’identifiant.

On peut maintenant représenter les données du prix médian de l’immobilier par EPCI sous forme de carte :

mf_map(x = data_immo, 
       var = "prix_med", 
       type = "choro",
       breaks = "quantile",
       nbreaks = 7,
       pal = "Mint",
       lwd = 0.01,
       leg_title = "Discrétisation par quantile",
       leg_val_rnd = 0)
mf_title("Prix médian de l'immobilier au m² par EPCI")
mf_credits("Sources: Notaires de France, IGN Admin Express")

Et avoir un aperçu rapide des autres données issues de l’INSEE :

par(mfrow = c(3,3))
for (var in colnames(data_immo)[6:13]) {
  mf_map(x = data_immo,
         var = var,
         type = "choro",
         breaks = "quantile",
         nbreaks = 7,
         pal = "Purples",
         lwd = 0.01,
         leg_pos = NA)
  mf_title(var)
  mf_credits('Sources: INSEE, IGN Admin Express')
}

Dans le cas où vous préféreriez manipuler vos données sous un format sp (package sp), ou dans le cas où ce format serait requis pour utiliser certains packages ou certaines formules, vous pouvez convertir votre objet sf en sp à l’aide de la ligne de commande suivante (nous en aurons besoin pour la suite) :

data_immo_sp <- as(data_immo, "Spatial")

Les packages sf (Simple Features for R) et sp (Classes and methods for spatial data) proposent tous deux des fonctions pour manipuler des données spatiales. Le package sf est plus récent et plus performant, mais beaucoup de packages R ne fonctionnent encore qu’avec le format sp.

2 Régression géographiquement pondérée (GWR)

Avant d’attaquer la GWR un petit appel des épisodes précédents !

Notre objectif est donc d’expliquer le prix median de l’immobilier des EPCI en France métropolitaine, en fonction d’un certian nombre d’autre variables comme par exemple la densité de population, le pourcentage de petits logement… Alors que nous pensions pouvoir réaliser une régression classique nous nous sommes rendu compte que cela ne fonctionnait pas ! Le spatial comme à son habitude venait perturbé nos plans.

Avant de poursuivre prenons le temps de visualiter la distribution de nos données le prix de l’immobilier et ses prédicteurs :

# Distribution de la variable dépendante :
add_histogram(plot_ly(data_immo, x = ~prix_med))

# Distribution des variables indépendantes :
a <- add_histogram(plot_ly(data_immo, x = ~log(perc_log_vac), name = "perc_log_vac"))
b <- add_histogram(plot_ly(data_immo, x = ~log(perc_maison), name = "perc_maison"))
c <- add_histogram(plot_ly(data_immo, x = ~log(perc_tiny_log), name = "perc_tiny_log"))
d <- add_histogram(plot_ly(data_immo, x = ~log(dens_pop), name = "dens_pop"))
e <- add_histogram(plot_ly(data_immo, x = ~log(med_niveau_vis), name = "med_niveau_vis"))
f <- add_histogram(plot_ly(data_immo, x = ~log(part_log_suroccup), name = "part_log_suroccup"))
g <- add_histogram(plot_ly(data_immo, x = ~log(part_agri_nb_emploi), name = "part_agri_nb_emploi"))
h <- add_histogram(plot_ly(data_immo, x = ~log(part_cadre_profintellec_nbemploi), name = "part_cadre_profintellec_nbemploi"))
fig = subplot(a, b, c, d, e, f, g, h, nrows = 2)
fig

Voici à quoi ressemble notre modèle linéaire classique

# Dans le fonctionnement sur R il est important de stocker la régression dans un objet.
# Pour lancer la régression on va utiliser la fonction lm() dont les 2 lettres sont l'acronyme pour linear model
mod.lm <- lm(formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + part_cadre_profintellec_nbemploi, 
             data = data_immo)

# On affiche les principaux résultats avec la fonction summary
summary(mod.lm)

## 
## Call:
## lm(formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + 
##     dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + 
##     part_cadre_profintellec_nbemploi, data = data_immo)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1566.8  -220.2   -27.2   174.0  3277.3 
## 
## Coefficients:
##                                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      1592.873     11.204 142.171  < 2e-16 ***
## perc_log_vac                     -287.337     14.134 -20.330  < 2e-16 ***
## perc_maison                      -128.255     20.041  -6.400 2.22e-10 ***
## perc_tiny_log                    -261.556     27.934  -9.363  < 2e-16 ***
## dens_pop                          173.942     15.272  11.389  < 2e-16 ***
## med_niveau_vis                    288.017     13.860  20.781  < 2e-16 ***
## part_log_suroccup                 388.982     27.352  14.221  < 2e-16 ***
## part_agri_nb_emploi               -20.785     15.059  -1.380    0.168    
## part_cadre_profintellec_nbemploi   -7.841     19.904  -0.394    0.694    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 391.8 on 1214 degrees of freedom
## Multiple R-squared:  0.7784, Adjusted R-squared:  0.7769 
## F-statistic: 532.9 on 8 and 1214 DF,  p-value: < 2.2e-16

mod.lm %>%
  tbl_regression(intercept = TRUE)

Characteristic	Beta	95% CI¹	p-value
(Intercept)	1,593	1,571, 1,615	<0.001
perc_log_vac	-287	-315, -260	<0.001
perc_maison	-128	-168, -89	<0.001
perc_tiny_log	-262	-316, -207	<0.001
dens_pop	174	144, 204	<0.001
med_niveau_vis	288	261, 315	<0.001
part_log_suroccup	389	335, 443	<0.001
part_agri_nb_emploi	-21	-50, 8.8	0.2
part_cadre_profintellec_nbemploi	-7.8	-47, 31	0.7
¹ CI = Confidence Interval

On peut également visualiser graphiquement les coefficients des variables explicatives avec le package GGally :

GGally::ggcoef_model(mod.lm)

Une manière très efficace de vérifier l’importance de la dimension spatiale c’est de cartographier nos résidus. S’il semble émerger une structure c’est que très vraissemblablement le spatial va jouer et que les conditions liées aux résidus de la réalisation d’une régression classique ne sont pas vérifiées.

De manière générale lorsque l’on travaille avec des données spatiales il est toujours bon de les cartographier.

On intègre les résidus à la table attributaire de notre objet sf. A priori, comme on a utilisé nos données spatiales (sf) pour la régression les données sont classées dans le bon ordre.

data_immo$res_reg <- mod.lm$residuals

On peut maintenant faire la carte des résidus :

# on fonce légèrement la couleur de fond pour mieux voir la classe centrale
mf_theme("default", bg = "#dedede")
# création d'une palette divergente avec une couleur neutre centrale
palette = mf_get_pal(n = c(nb_cl_inf0_res, nb_cl_sup0_res), pal = c("Teal", "Peach"), neutral = "#f5f5f5")
# la carte :
mf_map(x = data_immo, 
       var = "res_reg", 
       type = "choro", 
       border = NA,
       lwd = 0.1,
       breaks = breaks_residus,
       pal = palette,
       leg_title = "Valeur centrale =  0\nIntervalle = σ / 2", 
       leg_val_rnd = 1)
mf_title("Résidus de régression linéaire classique")
mf_credits("Sources: Notaires de France, INSEE, IGN Admin Express")

# réinitialisation du thème
mf_theme("default")

Sur cette carte on voit très clairement une spatialisation des résidus, sans même faire les tests nous aurions pu voir que la dimension spatiale jouait bien un rôle. Sans autocorrelation nous aurions eu une répartition aléatoire des résidus.

Comme nous l’avons vu, les données, lorsqu’elles sont spatialisées, sont souvent soumises aux phénomènes de dépendance et d’hétérogénéité spatiale.

Pour prendre en compte le phénomène de dépendance spatiale, on fait souvent appel aux régressions spatiales. Elles permettent à la fois de mieux comprendre la relation qui unit les variables explicatives à la variable étudiée d’une part ; et d’autre part de traiter la relation de notre VD à son propre voisinage. Les régressions spatiales peuvent donc (en fonction du type de régression choisie) intégrer les caractéristiques du voisinage (de la VD ou des VI) pour expliquer la VD. Il existe différents modèles de régression spatiale, toute la question est de savoir quel modèle utiliser ? Ce choix va dépendre de la nature des phénomènes spatiaux étudiés. Nous ne traiterons pas ici de ce sujet.

Par ailleurs, l’hétérogénéité, qui renvoie à une instabilité, induit une variabilité spatiale de nos paramètres. L’idée est que nos VI peuvent avoir un effet qui n’est pas le même partout dans l’espace. Dans ce cas nous optons pour la régression géographiquement pondérée (GWR).

Ainsi, si on s’intéresse au lien entre les voisins on est dans l’autocorrélation spatiale et les modèle SEM, SAR & cie mais si on s’intéresse à l’hétérogénéité de nos variables c’est-à-dire à leur variabilité selon leur localisation on est dans la GWR.

Attention, ce qui en théorie peut paraître assez tranché ne l’est souvent pas du tout en pratique. En effet, il y a bien des cas où l’on a du mal à savoir dans quel cadre on se situe exactement.

Pour réaliser une GWR sur R plusieurs packages reconnus existent. On peut citer notamment le package spgwr et le package GWmodel. Nous choisirons d’utiliser ici le package GWmodel.

2.1 Calcul de la matrice des distance

La première étape est de calculer la distance entre toutes nos observations. Pour ce faire nous utiliserons la fonction gw.dist().

# Le package GWmodel n'est pas compatible avec le format sf il a besoin d'un objet sp
# (contrairement à spgwr qui peut travailler avec un format sf)

# Pour construire la matrice de distances entre centroïdes des EPCI :
dm.calib <- gw.dist(dp.locat = coordinates(data_immo_sp))

2.2 Définition de la bande passante

La bande passante est une distance au-delà de laquelle le poids des observations est considéré comme nul. Le calcul de cette distance est très important car la valeur de la bande passante pourra fortement influencer notre modèle. La définition de la bande passante renvoie à quel type de pondération nous souhaitons appliquer. Heureusement la fonction bw.gwr va choisir pour nous le résultat optimal…

Pour ce faire la fonction va se baser sur un critère statistique que l’utilisateur devra définir : le CV (validation croisée) ou le AIC (Critère d’information d’Akaike). Elle reposera aussi sur un type de noyau qu’il faudra également définir : Gaussien, Exponentiel, Bicarré, Tricube ou encore Boxcar. Enfin, il sera également nécessaire de savoir si ce noyau pourra être adaptatif ou fixe.

Voici quelques informations pour guider nos choix :

Le critère CV a pour objectif de maximiser le pouvoir prédictif du modèle, le critère AIC va chercher un compromis entre le pouvoir prédictif du modèle et son degré de complexité. En général, le critère AIC est privilégié.
Avec un noyau fixe l’étendue du noyau est déterminée par la distance au point d’intérêt et il est identique en tout point de l’espace. Un noyau fixe est adapté si la répartition des données est homogène dans l’espace, l’unité de la bande passante sera donc une distance. Avec un noyau adaptatif l’étendue du noyau est déterminée par le nombre de voisins. Il est donc plus adapté à une répartition non homogène, l’unité sera alors le nombre de voisins.

Concernant la forme des noyaux :

Les noyaux gaussiens et exponentiels vont pondérer toutes les observations avec un poids qui tend vers zéro avec la distance au point estimé.
Les noyaux bisquare et tricube (dont les formes sont très proches) accordent également aux observations un poids décroissant avec la distance, mais par contre ce poids est nul au delà de la distance définie par la bande passante.
Le noyau Box-Car traite un phénomène continu de façon discontinue.

Manuel de géographie quantitative [@feuillet_2019]

Manuel de géographie quantitative (Feuillet et al. 2019)

Sachant que sur la forme du noyau, il est tout à fait possible de comparer deux pondérations et deux modèles de GWR.

Ici, nous allons tester avec un noyau gaussien, ce qui sera justifié un petit peu plus bas.

# Définition de la bande passante (bandwidth en anglais) :
bw_g <- bw.gwr(data = data_immo_sp, 
              approach = "AICc", 
              kernel = "gaussian", 
              adaptive = TRUE, 
              dMat = dm.calib,
              formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + part_cadre_profintellec_nbemploi)

## Adaptive bandwidth (number of nearest neighbours): 763 AICc value: 17959.04 
## Adaptive bandwidth (number of nearest neighbours): 480 AICc value: 17884.14 
## Adaptive bandwidth (number of nearest neighbours): 303 AICc value: 17798.07 
## Adaptive bandwidth (number of nearest neighbours): 196 AICc value: 17715.24 
## Adaptive bandwidth (number of nearest neighbours): 127 AICc value: 17622.76 
## Adaptive bandwidth (number of nearest neighbours): 87 AICc value: 17526.89 
## Adaptive bandwidth (number of nearest neighbours): 59 AICc value: 17422.28 
## Adaptive bandwidth (number of nearest neighbours): 45 AICc value: 17353.68 
## Adaptive bandwidth (number of nearest neighbours): 33 AICc value: 17283.09 
## Adaptive bandwidth (number of nearest neighbours): 29 AICc value: 17261.83 
## Adaptive bandwidth (number of nearest neighbours): 23 AICc value: 17250.23 
## Adaptive bandwidth (number of nearest neighbours): 22 AICc value: 17247.23 
## Adaptive bandwidth (number of nearest neighbours): 19 AICc value: 17201.79 
## Adaptive bandwidth (number of nearest neighbours): 19 AICc value: 17201.79

bw_g

## [1] 19

Notre bande passante est donc ici de 19 voisins, ce qui implique que les EPCI au-delà de cette “distance” auront un poids ramené à 0 et ne joueront donc plus de rôle dans la description de la relation statistique.

2.3 Estimation du modèle

Une fois la bande passante définie on peut lancer l’estimation de notre modèle de GWR :

mod.gwr_g <- gwr.robust(data = data_immo_sp, 
                   dMat = dm.calib,
                   bw = bw_g,
                   kernel = "gaussian",
                   filtered = FALSE, # un des problèmes de la GWR est de gérer des individus "aberrants" au niveau local. 2 méthodes ont été définies pour gérer cela. 
                                    # Méthode 1 (argument TRUE) on filtre en fonction des individus standardisés. L'objectif est de détecter les individus dont les résidus sont très élevés et de les exclure.
                                    # Méthode 2 (argument FALSE) on diminue le poids des observations aux résidus élevés.
                   adaptive = TRUE,
                   formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + part_cadre_profintellec_nbemploi)

Si on souhaite comparer deux modèles car nous avons un doute sur les paramètres c’est tout à fait possible. Par exemple ici nous souhaitons comparer deux formes de noyau :

bw_tri <- bw.gwr(data = data_immo_sp, 
              approach = "AICc", 
              kernel = "tricube", 
              adaptive = TRUE, 
              dMat = dm.calib,
              formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + part_cadre_profintellec_nbemploi)

## Adaptive bandwidth (number of nearest neighbours): 763 AICc value: 17701.17 
## Adaptive bandwidth (number of nearest neighbours): 480 AICc value: 17586.94 
## Adaptive bandwidth (number of nearest neighbours): 303 AICc value: 17422.51 
## Adaptive bandwidth (number of nearest neighbours): 196 AICc value: 17294.75 
## Adaptive bandwidth (number of nearest neighbours): 127 AICc value: 17254.1 
## Adaptive bandwidth (number of nearest neighbours): 87 AICc value: 17279.33 
## Adaptive bandwidth (number of nearest neighbours): 154 AICc value: 17259.05 
## Adaptive bandwidth (number of nearest neighbours): 112 AICc value: 17257.17 
## Adaptive bandwidth (number of nearest neighbours): 138 AICc value: 17254.8 
## Adaptive bandwidth (number of nearest neighbours): 122 AICc value: 17253.75 
## Adaptive bandwidth (number of nearest neighbours): 117 AICc value: 17255.04 
## Adaptive bandwidth (number of nearest neighbours): 123 AICc value: 17253.57 
## Adaptive bandwidth (number of nearest neighbours): 126 AICc value: 17254.25 
## Adaptive bandwidth (number of nearest neighbours): 123 AICc value: 17253.57

mod.gwr_tri <- gwr.robust(data = data_immo_sp, 
                   dMat = dm.calib,
                   bw = bw_tri,
                   kernel = "gaussian",
                   filtered = FALSE,
                   adaptive = TRUE,
                   formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + part_cadre_profintellec_nbemploi)

Best_gwr <- cbind(
  rbind(bw_g, bw_tri),
  rbind(mod.gwr_g$GW.diagnostic$gw.R2,mod.gwr_tri$GW.diagnostic$gw.R2),
  rbind(mod.gwr_g$GW.diagnostic$AIC,mod.gwr_tri$GW.diagnostic$AIC)) %>% 
  `colnames<-`(c("Nb Voisins","R2","AIC")) %>% 
  `rownames<-`(c("GAUSSIAN","TRICUBE"))

Best_gwr

##          Nb Voisins        R2      AIC
## GAUSSIAN         19 0.9324328 16849.26
## TRICUBE         123 0.8577370 17567.05

Le modèle avec une forme qui a été définie au format gaussien explique un meilleur \(R^2\) et le score d’\(AIC\) est plus faible. Ce modèle est donc plus performant et dans notre situation c’est plutôt ce modèle qu’il faut privilégier.

2.4 Interprétation des premiers résultats

Comme pour le modèle linéaire classique, l’objet qui contient notre GWR est composé de plusieurs éléments. Pour obtenir nos résultats il suffit d’appeler l’objet.

# Pour voir les différent élément qui compose notre modèle de GWR
summary(mod.gwr_g)

##               Length Class                    Mode
## GW.arguments    11   -none-                   list
## GW.diagnostic    8   -none-                   list
## lm              14   lm                       list
## SDF           1223   SpatialPolygonsDataFrame S4  
## timings          5   -none-                   list
## this.call       13   -none-                   call
## Ftests           0   -none-                   list

# Pour accéder aux résultat
mod.gwr_g

##    ***********************************************************************
##    *                       Package   GWmodel                             *
##    ***********************************************************************
##    Program starts at: 2023-10-07 14:25:49.656371 
##    Call:
##    gwr.basic(formula = formula, data = data, bw = bw, kernel = kernel, 
##     adaptive = adaptive, p = p, theta = theta, longlat = longlat, 
##     dMat = dMat, F123.test = F123.test, cv = T, W.vect = W.vect)
## 
##    Dependent (y) variable:  prix_med
##    Independent variables:  perc_log_vac perc_maison perc_tiny_log dens_pop med_niveau_vis part_log_suroccup part_agri_nb_emploi part_cadre_profintellec_nbemploi
##    Number of data points: 1223
##    ***********************************************************************
##    *                    Results of Global Regression                     *
##    ***********************************************************************
## 
##    Call:
##     lm(formula = formula, data = data)
## 
##    Residuals:
##     Min      1Q  Median      3Q     Max 
## -1566.8  -220.2   -27.2   174.0  3277.3 
## 
##    Coefficients:
##                                     Estimate Std. Error t value Pr(>|t|)    
##    (Intercept)                      1592.873     11.204 142.171  < 2e-16 ***
##    perc_log_vac                     -287.337     14.134 -20.330  < 2e-16 ***
##    perc_maison                      -128.255     20.041  -6.400 2.22e-10 ***
##    perc_tiny_log                    -261.556     27.934  -9.363  < 2e-16 ***
##    dens_pop                          173.942     15.272  11.389  < 2e-16 ***
##    med_niveau_vis                    288.017     13.860  20.781  < 2e-16 ***
##    part_log_suroccup                 388.982     27.352  14.221  < 2e-16 ***
##    part_agri_nb_emploi               -20.785     15.059  -1.380    0.168    
##    part_cadre_profintellec_nbemploi   -7.841     19.904  -0.394    0.694    
## 
##    ---Significance stars
##    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
##    Residual standard error: 391.8 on 1214 degrees of freedom
##    Multiple R-squared: 0.7784
##    Adjusted R-squared: 0.7769 
##    F-statistic: 532.9 on 8 and 1214 DF,  p-value: < 2.2e-16 
##    ***Extra Diagnostic information
##    Residual sum of squares: 186354789
##    Sigma(hat): 390.6721
##    AIC:  18086.13
##    AICc:  18086.31
##    BIC:  16985.31
##    ***********************************************************************
##    *          Results of Geographically Weighted Regression              *
##    ***********************************************************************
## 
##    *********************Model calibration information*********************
##    Kernel function: gaussian 
##    Adaptive bandwidth: 19 (number of nearest neighbours)
##    Regression points: the same locations as observations are used.
##    Distance metric: A distance matrix is specified for this model calibration.
## 
##    ****************Summary of GWR coefficient estimates:******************
##                                            Min.     1st Qu.      Median
##    Intercept                         1109.13394  1309.99018  1516.79508
##    perc_log_vac                      -919.82407  -224.60028  -163.40316
##    perc_maison                       -898.16365  -258.64481  -110.20592
##    perc_tiny_log                    -1210.96882  -206.68578   -92.91874
##    dens_pop                          -411.20172    72.82448   217.03593
##    med_niveau_vis                      81.29015   287.78992   358.32914
##    part_log_suroccup                 -543.68600    42.42839   142.48859
##    part_agri_nb_emploi               -498.74706   -65.81443   -32.88823
##    part_cadre_profintellec_nbemploi  -347.37935   -48.57329     0.58811
##                                         3rd Qu.    Max.
##    Intercept                         1696.07367 2330.78
##    perc_log_vac                      -113.83881  388.64
##    perc_maison                         76.01079  896.95
##    perc_tiny_log                       32.81544  773.30
##    dens_pop                           268.35651  659.84
##    med_niveau_vis                     413.36560  853.98
##    part_log_suroccup                  247.33650  700.50
##    part_agri_nb_emploi                 -1.15602  120.33
##    part_cadre_profintellec_nbemploi    49.93194  388.50
##    ************************Diagnostic information*************************
##    Number of data points: 1223 
##    Effective number of parameters (2trace(S) - trace(S'S)): 320.246 
##    Effective degrees of freedom (n-2trace(S) + trace(S'S)): 902.754 
##    AICc (GWR book, Fotheringham, et al. 2002, p. 61, eq 2.33): 17201.79 
##    AIC (GWR book, Fotheringham, et al. 2002,GWR p. 96, eq. 4.22): 16849.26 
##    BIC (GWR book, Fotheringham, et al. 2002,GWR p. 61, eq. 2.34): 17068.01 
##    Residual sum of squares: 56808914 
##    R-square value:  0.9324328 
##    Adjusted R-square value:  0.9084372 
## 
##    ***********************************************************************
##    Program stops at: 2023-10-07 14:26:25.464764

Cette visualisation des résultats nous propose d’abord un rappel complet du modèle linéaire classique. Puis viennent ensuite les informations concernant notre GWR. Le premier indicateur à analyser est le \(R^2 ajusté\) ajusté de la GWR, qui est nettement meilleur que celui de la régression linéaire multiple. On passe de \(77\%\) de variance expliqué à \(91\%\) avec la GWR. Ce \(R^2\), pour la GWR, est en fait une moyenne calculée des \(R^2\) de toutes les régressions locales réalisées par la GWR.

La seconde information qui nous intéresse particulièrement est les coefficient associés à nos VI. Nous voyons qu’ils ne sont pas présentés de la même manière que ceux de la régression linéaire. En effet, chaque VI va avoir des coefficients en fonction du minimum, maximum et des quartiles. Cela permet de rendre compte de la stationnarité de l’effet ou non. Dans notre cas on voit qu’il y a bien une variation et même dans certains cas une inversion des signes. Cela laisse supposer une non stationnarité des effets : un effet local peut être présent qui ne suivrait pas l’effet global.

Par exemple, pour le pourcentage de logement vacant avec un coefficient global (modèle linéaire) de \(-287\), quand ce pourcentage augmente le prix médian baisse. En simplifiant le pourcentage baisse d’une unité le prix médian augmente de \(287€\). Dans le cas de la densité de population on a un coefficient global positif donc une relation positive. La densité augmente donc le prix médian augmente. Ici, au global la densité augmente d’une unité le prix médian augmente de \(173€\).

Les résultats de la GWR peuvent donc être lus à une échelle globale pour mesurer la pertinence du modèle ; mais également à des échelles locales : les résultats illustrent ainsi comment les coefficients varient en fonction des unités spatiales. Gardons l’exemple de la densité de population. Dans les lieux où le prix médian est à son minimum le coefficient est de \(-411\) ; on a donc une relation négative. Dans ces espaces la densité augmente d’une unité le prix médian baisse de \(411€\).

Ensuite nous pouvons constater une inversion du signe du coefficient. Ainsi dans les EPCI du dernier quartile où le prix médian du logement est le plus élevé (par ex. Paris) le coefficient est positif. À son maximum une augmentation d’une unité entraîne une augmentation du prix de \(659€\). On a donc très clairement un effet de la densité qui ne sera pas du tout le même en fonction du lieu.

Nous pouvons également étudier l’intervalle interquartile. Ainsi, toujours pour la densité, ce résultat implique que pour 50% de nos unités spatiales (EPCI entre quartile 1 et 3), une augmentation d’une unité de la densité va impliquer une augmentation du prix médian entre \(72€\) et \(268€\).

Au travers de ces résultats on voit parfaitement comment une même variable peut avoir un effet différent, voire opposé en fonction des unités de lieu.

La cartographie va être la meilleure manière de représenter les betas (coefficients) et les différents indicateurs fournis avec la GWR, cela nous permet de décrire plus finement et plus précisément les phénomènes observés.

L’ensemble des données est stocké dans le sous objet SDF de notre modèle. Il contient l’ensemble des informations du modèle associé à chaque donnée spatiale.

On peut le convertir en un dataframe pour le visualiser plus facilement. À l’origine il est au format “SpatialPointsDataFrame”.

# Pour visualiser ce fichier dans R
#View(mod.gwr_g$SDF@data)

#Pour voir à quoi il ressemble dans ce document
datatable(mod.gwr_g$SDF@data)

# Pour voir les variables qui le constituent
names(mod.gwr_g$SDF@data)

##  [1] "Intercept"                           "perc_log_vac"                       
##  [3] "perc_maison"                         "perc_tiny_log"                      
##  [5] "dens_pop"                            "med_niveau_vis"                     
##  [7] "part_log_suroccup"                   "part_agri_nb_emploi"                
##  [9] "part_cadre_profintellec_nbemploi"    "y"                                  
## [11] "yhat"                                "residual"                           
## [13] "CV_Score"                            "Stud_residual"                      
## [15] "Intercept_SE"                        "perc_log_vac_SE"                    
## [17] "perc_maison_SE"                      "perc_tiny_log_SE"                   
## [19] "dens_pop_SE"                         "med_niveau_vis_SE"                  
## [21] "part_log_suroccup_SE"                "part_agri_nb_emploi_SE"             
## [23] "part_cadre_profintellec_nbemploi_SE" "Intercept_TV"                       
## [25] "perc_log_vac_TV"                     "perc_maison_TV"                     
## [27] "perc_tiny_log_TV"                    "dens_pop_TV"                        
## [29] "med_niveau_vis_TV"                   "part_log_suroccup_TV"               
## [31] "part_agri_nb_emploi_TV"              "part_cadre_profintellec_nbemploi_TV"
## [33] "E_weigts"                            "Local_R2"

# Intercept : c'est la constante c'est à dire prix médian de référence
# nom de la variable : estimation du coefficient, du beta associé à la VI en chaque point.
# y : les valeurs de la VD
# yhat : valeur de y prédite.
# residual, Stud_residual : résidu et résidu standardisé
# CV_score : score de validation croisée
# _SE : erreur standard de l'estimation du coefficient
# _TV : t-value de l'estimation du coefficient
# E_weight : poids des observations dans la régression robuste
# Local_R2 : R2 au niveau de chaque unité spatiale

2.4.1 Étude des résidus

Commençons par une étude des résidus afin de vérifier que cette fois ils n’ont pas de structure apparente.

# on récupère les résidus dans data_immo
res_gwr <- mod.gwr_g$SDF$Stud_residual
data_immo$res_gwr <- res_gwr
# calcul des limites de classes avec la fonction discr, centrées sur 0
res_resgwr <- discr(data_immo$res_gwr, 0, "class_center", sd(data_immo$res_gwr)*0.5, 10)
breaks_gwr <- res_resgwr[[1]]
nb_cl_sup0_gwr <- res_resgwr[[2]]
nb_cl_inf0_gwr <- res_resgwr[[3]]
# création de la palette correspondante
palette = mf_get_pal(n = c(nb_cl_inf0_gwr, nb_cl_sup0_gwr), pal = c("Teal", "Peach"), neutral = "#f5f5f5")
# la carte des résidus
mf_map(x = data_immo, 
       var = "res_gwr", 
       type = "choro", 
       border = "gray", 
       lwd = 0.2, 
       breaks = breaks_gwr,
       pal = palette, 
       leg_title = "Discrétisation standardisée :\nvaleur centrale = 0\nintervalle = σ / 2", 
       leg_val_rnd = 1)
mf_title("Résidus GWR")
mf_credits("Sources: Notaires de France, INSEE, IGN Admin Express")

Cette carte ne présente pas de structure spatiale marquée et nous amène à penser que nous avons expliqué l’ensemble des phénomènes spatiaux liés aux questions de prix de l’immobilier.

2.4.2 Étude des coefficients

Pour visualiser la non stationnarité des effets de nos VI la solution la plus efficace est la carte.

# On ajoute à data_immo les coefficients
data_immo$agri.coef <- mod.gwr_g$SDF$part_agri_nb_emploi
data_immo$perc_maison.coef <- mod.gwr_g$SDF$perc_maison
data_immo$dens_pop.coef <- mod.gwr_g$SDF$dens_pop
data_immo$med_vie.coef <- mod.gwr_g$SDF$med_niveau_vis
data_immo$logvac.coef <- mod.gwr_g$SDF$perc_log_vac
data_immo$tinylog.coef <- mod.gwr_g$SDF$perc_tiny_log
data_immo$suroccup.coef <- mod.gwr_g$SDF$part_log_suroccup
data_immo$cadre.coef <- mod.gwr_g$SDF$part_cadre_profintellec_nbemploi

Pour réaliser les cartes, avec une discrétisation standardisée centrée sur 0 :

par(mfrow = c(4, 2)) 
for (var in colnames(data_immo)[17:24]) {
  # calcul des limites de classe
  res <- discr(data.frame(data_immo)[, var], 0, "class_center", sd(data.frame(data_immo)[, var])*0.5, 10)
  breaks <- res[[1]]
  # palette de couleurs
  nb_cl_sup0 <- res[[2]]
  nb_cl_inf0 <- res[[3]]
  if (nb_cl_inf0 > 0) {
    palette = mf_get_pal(n = c(nb_cl_inf0, nb_cl_sup0), pal = c("Teal", "Peach"), neutral = "#f5f5f5")
  } else { # cas de la médiane du niveau de vie où la valeur min est supérieure à 0
    palette = mf_get_pal(n = c(nb_cl_sup0), pal = c("Peach"))
  }
  # la carte
  mf_map(x = data_immo,
         var = var,
         type = "choro",
         border = "gray",
         lwd = 0.1,
         breaks = breaks,
         pal = palette,
         leg_pos = "left",
         leg_title = NA,
         leg_val_rnd = 0)
  mf_title(var)
}

Les cartes des betas vont illustrer la variation des effets en fonctions des entités spatiales et de leur voisinage. Dans notre cas on verra quels sont les EPCI où l’effet du coefficient est négatif et ceux où il est positif, c’est-à-dire dans quel EPCI notre VI va entraîner une augmentation du prix médian et dans quel autre au contraire une diminution, toutes choses égales par ailleurs. Sachant que dans notre cas toutes les VI sont significatives, elles ont donc toutes un effet qui varie localement.

Au-delà de cette visualisation VI par VI, il peut être intéressant de savoir par EPCI quelle variable sera la plus explicative dans la relation à notre VD, laquelle a l’impact le plus important. Nous avons donc réalisé une carte des contributions max par EPCI. Pour la réaliser nous nous sommes basés sur le T-value.

data_immo$agri.t <- mod.gwr_g$SDF$part_agri_nb_emploi_TV
data_immo$maison.t <- mod.gwr_g$SDF$perc_maison_TV
data_immo$dens.t <- mod.gwr_g$SDF$dens_pop_TV
data_immo$medvie.t <- mod.gwr_g$SDF$med_niveau_vis_TV
data_immo$logvac.t <- mod.gwr_g$SDF$perc_log_vac_TV
data_immo$tinylog.t <- mod.gwr_g$SDF$perc_tiny_log_TV
data_immo$suroccup.t <- mod.gwr_g$SDF$part_log_suroccup_TV
data_immo$cadre.t <- mod.gwr_g$SDF$part_cadre_profintellec_nbemploi_TV     

# Définir contrib max
df <- as.data.frame(data_immo)
# On passe les t-values en valeurs absolues pour voir la plus grande contribution dans un sens sens ou dans l'autre
data_immo$contribmax<- colnames(df[, c(25:32)])[max.col(abs(df[, c(25:32)]),ties.method="first")]

par(mfrow = c(1, 1))
# Carte
mf_map(x = data_immo, 
       var = "contribmax", 
       type = "typo", 
       pal = brewer.pal(6,'Set2'),
       border = "white",
       lwd = 0.2)
mf_title("Carte des variables contribuant le plus par epci")

Au-delà de la non-stationnarité de nos VI, cette carte met en évidence un autre phénomène : Dans un modèle linéaire classique, la sélection des VI se fait de manière itérative (ascendante, descendante, ou mixte). L’inclusion d’une VI a pour conséquence d’en exclure d’autres qui présenteraient trop de multi-colinéarité. Le phénomène mis en évidence ici repose sur le fait qu’en fonction du lieu, la hiérarchie des VI, et donc leur pertinence au sein du modèle, n’est pas constante. Cette visualisation des données ouvre donc une perspective intéressante pour la suite. Il serait tout à fait pertinent de développer une méthode permettant d’effectuer un stepwise préalable à toute GWR valable pour chaque lieu individuellement. Un tel modèle traiterait dans sa totalité le problème de non-stationnarité.

Nous pouvons également cartographier les \(R^2\) locaux, ce qui fournit une indication sur les zones où la variabilité est la mieux expliquée.

data_immo$r2local=mod.gwr_g$SDF$Local_R2

mf_map(x = data_immo, 
       var = "r2local", 
       type = "choro",
       breaks = "quantile",
       nbreaks = 11,
       pal= "Reds",
       border = "gray",
       lwd = 0.2,
       leg_title = "Discrétisation par quantile")
mf_title("R² locaux")

Nous avons choisi d’utiliser une palette de valeurs large pour représenter au mieux les différences de \(R^2\) locaux. Toutefois, cette carte peut paraître trompeuse ; tous les EPCI obtiennent une explication satisfaisante avec un R² local minimum de 0.68.

À partir des t-value on peut aussi étudier la significativité des effets sur le territoire. On peut ainsi calculer et cartographier un indicateur qui représenterait le nombre de VI dont l’effet est significatif sur chaque unité spatiale. Cela donne une bonne idée de la complexité du phénomène sur un espace donné (en effet sur un EPCI on peut avoir toutes les variables significatives, elle jouent donc sur cet espace toutes un rôles) et souligne l’importance d’avoir une carte par coefficient. Cette carte montre également qu’un modèle parfaitement adaptatif gagnerait en sobriété et donnerait un AIC plus satisfaisant en ne sélectionnant localement que les variables réellement significatives dans la relation.

# Pour rappel si on a plus de 200 individus et le t-value > |1.96| on pourra considérer le coefficient comme significatif au seuil de 0.05 (95% chances que ce ne soit pas dû au hasard)

data_immo$nbsignif_t <- rowSums(abs(df[, c(25:32)]) > 1.96)

mf_map(x = data_immo, 
       var = "nbsignif_t", 
       type = "typo",
       pal = "Reds",
       border = "gray",
       lwd = 0.2)
mf_title("Nombre de Betas significatifs par EPCI (t-value)")

Il se peut que cela soit plus intéressant d’utiliser les p-value, notamment si vous avez moins de 200 individus.

# Les p-value ne sont pas fournis dans le modèle de la GWR, on pourrait les calculer à partir de t-value et de l'erreur standard mais le package GWmodel propose une fonction pour les obtenir
pvalue <- gwr.t.adjust(mod.gwr_g)

# On ajoute les p-value à notre fichier
data_immo$agri.p <- pvalue$SDF$part_agri_nb_emploi_p 
data_immo$maison.p <- pvalue$SDF$perc_maison_p
data_immo$dens.p <- pvalue$SDF$dens_pop_p
data_immo$medvie.p <- pvalue$SDF$med_niveau_vis_p
data_immo$logvac.p <- pvalue$SDF$perc_log_vac_p
data_immo$tinylog.p <- pvalue$SDF$perc_tiny_log_p
data_immo$suroccup.p <- pvalue$SDF$part_log_suroccup_p
data_immo$cadre.p <- pvalue$SDF$part_cadre_profintellec_nbemploi_p

df<- as.data.frame(data_immo)
data_immo$nbsignif_p <- rowSums(df[, c(36:43)] < 0.05)

mf_map(x = data_immo, 
       var = "nbsignif_p", 
       type = "typo",
       pal= "Reds",
       border = "gray",
       lwd = 0.2,)
mf_title("Nombre des Betas significatifs par EPCI (p-value)")

Sur ces deux dernières cartes (nombre de betas significatifs par EPCI avec les t-value et p-value) nous avons pris volontairement une liberté avec les règles de sémiologie graphique de Bertin. En effet, s’agissant de valeurs discrètes, nous aurions dû les représenter sous forme de symboles proportionnels. Cependant, dans la mesure où les EPCI ont une taille relativement homogène et par souci de visibilité, et aussi parce qu’il s’agit plus d’une carte exploratoire que d’un véritable rendu, nous avons opté pour une carte choroplèthe. Que les mânes de Bertin nous pardonnent !

Dans ce cadre, il est possible de réaliser une collection de cartes des p-value (ou t-value) comme ce qui a été fait pour les coefficients. L’intérêt est de voir où l’effet de la VI est significatif et où il ne l’est pas.

# Ici nous représenterons les p-value avec un découpage par classe de significativité et seulement les p-value de 2 VI

par(mfrow = c(1, 2))

# Par exemple les p-value des coefficients de la variable part de l'emploi agriculteur
data_immo<- data_immo %>%  mutate(agri.p_fac = case_when(agri.p<= 0.002 ~ "[0;0.002[",
                           agri.p <= 0.01 ~ "[0.002;0.01[",
                           agri.p <= 0.05 ~ "[0.01;0.05[",
                           agri.p <= 0.1 ~ "[0.05;0.1[",
                           TRUE ~ "[0.1;1]")) %>%
  mutate(agri.p_fac = factor(agri.p_fac,
                        levels = c("[0;0.002[", "[0.002;0.01[",
                                   "[0.01;0.05[",
                                  "[0.05;0.1[", "[0.1;1]")))


mypal2 <- mf_get_pal(n = 5, palette = "OrRd")

mf_map(x = data_immo, 
       var = "agri.p_fac", 
       type = "typo", 
       border = "grey3", 
       lwd = 0.1, 
       pal = mypal2, 
       leg_title = "Classe P-value")
mf_title("P-value du coefficient de la part d'emploi agriculteurs")

# Pour la densité de population
data_immo<- data_immo %>%  mutate(dens.p_fac = case_when(dens.p <= 0.002 ~ "[0;0.002[",
                           dens.p <= 0.01 ~ "[0.002;0.01[",
                           dens.p <= 0.05 ~ "[0.01;0.05[",
                           dens.p <= 0.1 ~ "[0.05;0.1[",
                           TRUE ~ "[0.1;1]")) %>%
  mutate(dens.p_fac = factor(dens.p_fac,
                        levels = c("[0;0.002[", "[0.002;0.01[",
                                   "[0.01;0.05[",
                                  "[0.05;0.1[", "[0.1;1]")))

mypal2 <- mf_get_pal(n = 5, palette = "OrRd")

mf_map(x = data_immo, 
       var = "dens.p_fac", 
       type = "typo", 
       border = "grey3", 
       lwd = 0.1, 
       pal=mypal2, 
       leg_title = "Classe P-value")
mf_title("P-value du coefficient de la densité de population")

Ces cartes des p-value sont particulièrement importantes car elles nous donnent les endroits où l’effet est significatif. En effet, on sait que la VI a effet global qui est significatif, qu’elle a en plus une variabilité locale or localement elle n’est pas partout significative. Pour la part d’agriculteur dans l’emploi, l’effet est significatif quasiment uniquement dans le Sud-Est.

3 Pour aller plus loin…

Cette partie et le code qui va avec est directement tiré du cours donnée à SIGR par Thierry Feuillet (feuillet?).

3.1 GWR Multiscalaire

Selon Thierry Feuillet concernant la GWR multiscalaire :

“Il n’y a pas de raison de penser que tous les prédicteurs agissent sur le prix à la même échelle (c’est-à-dire selon un même schéma de voisinage). Certains processus peuvent être locaux, d’autres globaux. Récemment, une extension de la GWR a été proposée, permettant de relâcher cette hypothèse d’égalité des échelles : la GWR multiscalaire (MGWR, Fotheringham et al., 2017). Le principe est simple : un algorithme optimise le choix de la bandwidth pour chaque prédicteur, en fonction des autres. Il en résulte un modèle souvent mixte.”

Avec notre modèle même simplifié avec 3 prédicteurs le modèle est très long à tourner à l’échelle de la France entière, plus de 2h !

On part donc ici sur un corpus plus petit : La Bretagne sans toucher à nos VI et au nombre d’itérations

La GWR multiscalaire est une méthode itérative, la faire tourner peut donc être TRES long. Plusieurs possibilités : diminuer le nombre maximum d’itérations, simplifier le modèle et intégrer moins de VI ou sinon réduire l’emprise spatiale.

Afin de réduire les temps de traitement, on filtre d’abord les données pour ne garder que les EPCI bretons. Un EPCI pouvant être à cheval sur 2 régions, on ne va garder ici que les EPCI dont le centroïde est en région Bretagne. On va pour ça se servir de la couche REGION.shp de la base ADMIN-EXPRESS de l’IGN.

# Lecture de la couche région dans un objet sf
shp_path <- here("data", "REGION.shp")
region_sf <- st_read(shp_path)

## Reading layer `REGION' from data source 
##   `/media/glecampion/Data/gregoire/atelierR_foss_gwr/data/REGION.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 13 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 99040 ymin: 6046546 xmax: 1242445 ymax: 7110479
## Projected CRS: RGF93 v1 / Lambert-93

# Pour ne garder que la Bretagne
bzh_sf <- region_sf[region_sf$NOM == "Bretagne",]
mf_map(x = bzh_sf)

# Création des centroïdes des EPCI
epci_centroids <- st_centroid(epci_sf)

## Warning: st_centroid assumes attributes are constant over geometries

# Sélection des centroïdes dans la région Bretagne
epci_centroids_bzh <- st_intersection(epci_centroids, bzh_sf)

## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

mf_map(x = epci_centroids_bzh)

# On ne garde que les EPCI correspondant à ces centroïdes
epci_bzh <- merge(x = epci_sf, y = st_drop_geometry(epci_centroids_bzh), by.x = "CODE_SIREN", by.y = "CODE_SIREN")
mf_map(x = epci_bzh)

# Il faut refaire la jointure entre les données sur les EPCI et l'objet sf
data_immo_bzh <- merge(x = epci_bzh, y = immo_df, by.x = "CODE_SIREN", by.y = "SIREN")
# conversion objet sf vers objet sp pour le package GWmodel
data_immo_bzh_sp <-as(data_immo_bzh, "Spatial")

On peut maintenant lancer la GWR multiscalaire :

#source("gwr.multiscale_T.r")
# On lance la MGWR
# On note qu'il n'est pas ici nécessaire de définir une bande passante. Le principe de la GWR multiscalaire est justement d'adapter à chaque relation locale une bande passante (d'où les itérations)
MGWR <- gwr.multiscale(formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + part_cadre_profintellec_nbemploi, 
                       data = data_immo_bzh_sp, 
                       kernel = "gaussian", 
                       predictor.centered = rep(T, 3), # centrage des prédicteurs
                       adaptive = TRUE,
                       bws0 = rep(1,4)) # BW minimum pour l'optimisation

## Warning in gwr.multiscale(formula = prix_med ~ perc_log_vac + perc_maison + :
## All the predictors will be centered, please check the parameter
## predictor.centered

## ------   Calculate the initial beta0 from the above bandwidths    ------
## ------            The end for calculating the initial beta0              ------
## ------ Select the optimum bandwidths for each independent variable via  AIC  aproach ------
## *****  The back-fitting process for model calibration with bandwiths selected *****
##     Iteration  1 :
## Now select an optimum bandwidth for the variable:  Intercept 
## The newly selected bandwidth for variable  Intercept is:  14 
## The bandwidth used in the last ieration is:  1 and the difference between these two bandwidths is:  13 
## The bandwidth for variable  Intercept will be continually selected in the next ieration.
## Now select an optimum bandwidth for the variable:  perc_log_vac 
## The newly selected bandwidth for variable  perc_log_vac is:  53 
## The bandwidth used in the last ieration is:  1 and the difference between these two bandwidths is:  52 
## The bandwidth for variable  perc_log_vac will be continually selected in the next ieration.
## Now select an optimum bandwidth for the variable:  perc_maison 
## The newly selected bandwidth for variable  perc_maison is:  39 
## The bandwidth used in the last ieration is:  1 and the difference between these two bandwidths is:  38 
## The bandwidth for variable  perc_maison will be continually selected in the next ieration.
## Now select an optimum bandwidth for the variable:  perc_tiny_log 
## The newly selected bandwidth for variable  perc_tiny_log is:  39 
## The bandwidth used in the last ieration is:  1 and the difference between these two bandwidths is:  38 
## The bandwidth for variable  perc_tiny_log will be continually selected in the next ieration.
## Now select an optimum bandwidth for the variable:  dens_pop 
## The newly selected bandwidth for variable  dens_pop is:  54 
## The bandwidth used in the last ieration is:  1 and the difference between these two bandwidths is:  53 
## The bandwidth for variable  dens_pop will be continually selected in the next ieration.
## Now select an optimum bandwidth for the variable:  med_niveau_vis 
## The newly selected bandwidth for variable  med_niveau_vis is:  56 
## The bandwidth used in the last ieration is:  1 and the difference between these two bandwidths is:  55 
## The bandwidth for variable  med_niveau_vis will be continually selected in the next ieration.
## Now select an optimum bandwidth for the variable:  part_log_suroccup 
## The newly selected bandwidth for variable  part_log_suroccup is:  56 
## The bandwidth used in the last ieration is:  1 and the difference between these two bandwidths is:  55 
## The bandwidth for variable  part_log_suroccup will be continually selected in the next ieration.
## Now select an optimum bandwidth for the variable:  part_agri_nb_emploi 
## The newly selected bandwidth for variable  part_agri_nb_emploi is:  56 
## The bandwidth used in the last ieration is:  1 and the difference between these two bandwidths is:  55 
## The bandwidth for variable  part_agri_nb_emploi will be continually selected in the next ieration.
## Now select an optimum bandwidth for the variable:  part_cadre_profintellec_nbemploi 
## The newly selected bandwidth for variable  part_cadre_profintellec_nbemploi is:  56 
## The bandwidth used in the last ieration is:  1 and the difference between these two bandwidths is:  55 
## The bandwidth for variable  part_cadre_profintellec_nbemploi will be continually selected in the next ieration.
##     Ieration  1 the differential change value of RSS (dCVR) is:  0.4470374 
## ----------End of    Iteration  1 ----------
##     Iteration  2 :
## Now select an optimum bandwidth for the variable:  Intercept 
## The newly selected bandwidth for variable  Intercept is:  11 
## The bandwidth used in the last ieration is:  14 and the difference between these two bandwidths is:  3 
## The bandwidth for variable  Intercept will be continually selected in the next ieration.
## Now select an optimum bandwidth for the variable:  perc_log_vac 
## The newly selected bandwidth for variable  perc_log_vac is:  53 
## The bandwidth used in the last ieration is:  53 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  perc_log_vac seems to be converged for  1  times.It will be continually optimized in the next  4  times
## Now select an optimum bandwidth for the variable:  perc_maison 
## The newly selected bandwidth for variable  perc_maison is:  56 
## The bandwidth used in the last ieration is:  39 and the difference between these two bandwidths is:  17 
## The bandwidth for variable  perc_maison will be continually selected in the next ieration.
## Now select an optimum bandwidth for the variable:  perc_tiny_log 
## The newly selected bandwidth for variable  perc_tiny_log is:  39 
## The bandwidth used in the last ieration is:  39 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  perc_tiny_log seems to be converged for  1  times.It will be continually optimized in the next  4  times
## Now select an optimum bandwidth for the variable:  dens_pop 
## The newly selected bandwidth for variable  dens_pop is:  56 
## The bandwidth used in the last ieration is:  54 and the difference between these two bandwidths is:  2 
## The bandwidth for variable  dens_pop will be continually selected in the next ieration.
## Now select an optimum bandwidth for the variable:  med_niveau_vis 
## The newly selected bandwidth for variable  med_niveau_vis is:  56 
## The bandwidth used in the last ieration is:  56 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  med_niveau_vis seems to be converged for  1  times.It will be continually optimized in the next  4  times
## Now select an optimum bandwidth for the variable:  part_log_suroccup 
## The newly selected bandwidth for variable  part_log_suroccup is:  56 
## The bandwidth used in the last ieration is:  56 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  part_log_suroccup seems to be converged for  1  times.It will be continually optimized in the next  4  times
## Now select an optimum bandwidth for the variable:  part_agri_nb_emploi 
## The newly selected bandwidth for variable  part_agri_nb_emploi is:  56 
## The bandwidth used in the last ieration is:  56 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  part_agri_nb_emploi seems to be converged for  1  times.It will be continually optimized in the next  4  times
## Now select an optimum bandwidth for the variable:  part_cadre_profintellec_nbemploi 
## The newly selected bandwidth for variable  part_cadre_profintellec_nbemploi is:  56 
## The bandwidth used in the last ieration is:  56 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  part_cadre_profintellec_nbemploi seems to be converged for  1  times.It will be continually optimized in the next  4  times
##     Ieration  2 the differential change value of RSS (dCVR) is:  0.2210537 
## ----------End of    Iteration  2 ----------
##     Iteration  3 :
## Now select an optimum bandwidth for the variable:  Intercept 
## The newly selected bandwidth for variable  Intercept is:  11 
## The bandwidth used in the last ieration is:  11 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  Intercept seems to be converged for  1  times.It will be continually optimized in the next  4  times
## Now select an optimum bandwidth for the variable:  perc_log_vac 
## The newly selected bandwidth for variable  perc_log_vac is:  53 
## The bandwidth used in the last ieration is:  53 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  perc_log_vac seems to be converged for  2  times.It will be continually optimized in the next  3  times
## Now select an optimum bandwidth for the variable:  perc_maison 
## The newly selected bandwidth for variable  perc_maison is:  56 
## The bandwidth used in the last ieration is:  56 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  perc_maison seems to be converged for  1  times.It will be continually optimized in the next  4  times
## Now select an optimum bandwidth for the variable:  perc_tiny_log 
## The newly selected bandwidth for variable  perc_tiny_log is:  39 
## The bandwidth used in the last ieration is:  39 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  perc_tiny_log seems to be converged for  2  times.It will be continually optimized in the next  3  times
## Now select an optimum bandwidth for the variable:  dens_pop 
## The newly selected bandwidth for variable  dens_pop is:  56 
## The bandwidth used in the last ieration is:  56 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  dens_pop seems to be converged for  1  times.It will be continually optimized in the next  4  times
## Now select an optimum bandwidth for the variable:  med_niveau_vis 
## The newly selected bandwidth for variable  med_niveau_vis is:  56 
## The bandwidth used in the last ieration is:  56 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  med_niveau_vis seems to be converged for  2  times.It will be continually optimized in the next  3  times
## Now select an optimum bandwidth for the variable:  part_log_suroccup 
## The newly selected bandwidth for variable  part_log_suroccup is:  56 
## The bandwidth used in the last ieration is:  56 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  part_log_suroccup seems to be converged for  2  times.It will be continually optimized in the next  3  times
## Now select an optimum bandwidth for the variable:  part_agri_nb_emploi 
## The newly selected bandwidth for variable  part_agri_nb_emploi is:  56 
## The bandwidth used in the last ieration is:  56 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  part_agri_nb_emploi seems to be converged for  2  times.It will be continually optimized in the next  3  times
## Now select an optimum bandwidth for the variable:  part_cadre_profintellec_nbemploi 
## The newly selected bandwidth for variable  part_cadre_profintellec_nbemploi is:  56 
## The bandwidth used in the last ieration is:  56 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  part_cadre_profintellec_nbemploi seems to be converged for  2  times.It will be continually optimized in the next  3  times
##     Ieration  3 the differential change value of RSS (dCVR) is:  0.1316818 
## ----------End of    Iteration  3 ----------
##     Iteration  4 :
## Now select an optimum bandwidth for the variable:  Intercept 
## The newly selected bandwidth for variable  Intercept is:  11 
## The bandwidth used in the last ieration is:  11 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  Intercept seems to be converged for  2  times.It will be continually optimized in the next  3  times
## Now select an optimum bandwidth for the variable:  perc_log_vac 
## The newly selected bandwidth for variable  perc_log_vac is:  53 
## The bandwidth used in the last ieration is:  53 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  perc_log_vac seems to be converged for  3  times.It will be continually optimized in the next  2  times
## Now select an optimum bandwidth for the variable:  perc_maison 
## The newly selected bandwidth for variable  perc_maison is:  56 
## The bandwidth used in the last ieration is:  56 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  perc_maison seems to be converged for  2  times.It will be continually optimized in the next  3  times
## Now select an optimum bandwidth for the variable:  perc_tiny_log 
## The newly selected bandwidth for variable  perc_tiny_log is:  39 
## The bandwidth used in the last ieration is:  39 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  perc_tiny_log seems to be converged for  3  times.It will be continually optimized in the next  2  times
## Now select an optimum bandwidth for the variable:  dens_pop 
## The newly selected bandwidth for variable  dens_pop is:  54 
## The bandwidth used in the last ieration is:  56 and the difference between these two bandwidths is:  2 
## The bandwidth for variable  dens_pop will be continually selected in the next ieration.
## Now select an optimum bandwidth for the variable:  med_niveau_vis 
## The newly selected bandwidth for variable  med_niveau_vis is:  56 
## The bandwidth used in the last ieration is:  56 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  med_niveau_vis seems to be converged for  3  times.It will be continually optimized in the next  2  times
## Now select an optimum bandwidth for the variable:  part_log_suroccup 
## The newly selected bandwidth for variable  part_log_suroccup is:  56 
## The bandwidth used in the last ieration is:  56 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  part_log_suroccup seems to be converged for  3  times.It will be continually optimized in the next  2  times
## Now select an optimum bandwidth for the variable:  part_agri_nb_emploi 
## The newly selected bandwidth for variable  part_agri_nb_emploi is:  56 
## The bandwidth used in the last ieration is:  56 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  part_agri_nb_emploi seems to be converged for  3  times.It will be continually optimized in the next  2  times
## Now select an optimum bandwidth for the variable:  part_cadre_profintellec_nbemploi 
## The newly selected bandwidth for variable  part_cadre_profintellec_nbemploi is:  56 
## The bandwidth used in the last ieration is:  56 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  part_cadre_profintellec_nbemploi seems to be converged for  3  times.It will be continually optimized in the next  2  times
##     Ieration  4 the differential change value of RSS (dCVR) is:  0.1174028 
## ----------End of    Iteration  4 ----------
##     Iteration  5 :
## Now select an optimum bandwidth for the variable:  Intercept 
## The newly selected bandwidth for variable  Intercept is:  11 
## The bandwidth used in the last ieration is:  11 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  Intercept seems to be converged for  3  times.It will be continually optimized in the next  2  times
## Now select an optimum bandwidth for the variable:  perc_log_vac 
## The newly selected bandwidth for variable  perc_log_vac is:  53 
## The bandwidth used in the last ieration is:  53 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  perc_log_vac seems to be converged for  4  times.It will be continually optimized in the next  1  times
## Now select an optimum bandwidth for the variable:  perc_maison 
## The newly selected bandwidth for variable  perc_maison is:  56 
## The bandwidth used in the last ieration is:  56 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  perc_maison seems to be converged for  3  times.It will be continually optimized in the next  2  times
## Now select an optimum bandwidth for the variable:  perc_tiny_log 
## The newly selected bandwidth for variable  perc_tiny_log is:  39 
## The bandwidth used in the last ieration is:  39 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  perc_tiny_log seems to be converged for  4  times.It will be continually optimized in the next  1  times
## Now select an optimum bandwidth for the variable:  dens_pop 
## The newly selected bandwidth for variable  dens_pop is:  47 
## The bandwidth used in the last ieration is:  54 and the difference between these two bandwidths is:  7 
## The bandwidth for variable  dens_pop will be continually selected in the next ieration.
## Now select an optimum bandwidth for the variable:  med_niveau_vis 
## The newly selected bandwidth for variable  med_niveau_vis is:  56 
## The bandwidth used in the last ieration is:  56 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  med_niveau_vis seems to be converged for  4  times.It will be continually optimized in the next  1  times
## Now select an optimum bandwidth for the variable:  part_log_suroccup 
## The newly selected bandwidth for variable  part_log_suroccup is:  56 
## The bandwidth used in the last ieration is:  56 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  part_log_suroccup seems to be converged for  4  times.It will be continually optimized in the next  1  times
## Now select an optimum bandwidth for the variable:  part_agri_nb_emploi 
## The newly selected bandwidth for variable  part_agri_nb_emploi is:  56 
## The bandwidth used in the last ieration is:  56 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  part_agri_nb_emploi seems to be converged for  4  times.It will be continually optimized in the next  1  times
## Now select an optimum bandwidth for the variable:  part_cadre_profintellec_nbemploi 
## The newly selected bandwidth for variable  part_cadre_profintellec_nbemploi is:  56 
## The bandwidth used in the last ieration is:  56 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  part_cadre_profintellec_nbemploi seems to be converged for  4  times.It will be continually optimized in the next  1  times
##     Ieration  5 the differential change value of RSS (dCVR) is:  0.1086761 
## ----------End of    Iteration  5 ----------
##     Iteration  6 :
## Now select an optimum bandwidth for the variable:  Intercept 
## The newly selected bandwidth for variable  Intercept is:  11 
## The bandwidth used in the last ieration is:  11 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  Intercept seems to be converged for  4  times.It will be continually optimized in the next  1  times
## Now select an optimum bandwidth for the variable:  perc_log_vac 
## The newly selected bandwidth for variable  perc_log_vac is:  53 
## The bandwidth used in the last ieration is:  53 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  perc_log_vac seems to be converged and will be kept the same in the following ierations.
## Now select an optimum bandwidth for the variable:  perc_maison 
## The newly selected bandwidth for variable  perc_maison is:  56 
## The bandwidth used in the last ieration is:  56 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  perc_maison seems to be converged for  4  times.It will be continually optimized in the next  1  times
## Now select an optimum bandwidth for the variable:  perc_tiny_log 
## The newly selected bandwidth for variable  perc_tiny_log is:  39 
## The bandwidth used in the last ieration is:  39 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  perc_tiny_log seems to be converged and will be kept the same in the following ierations.
## Now select an optimum bandwidth for the variable:  dens_pop 
## The newly selected bandwidth for variable  dens_pop is:  41 
## The bandwidth used in the last ieration is:  47 and the difference between these two bandwidths is:  6 
## The bandwidth for variable  dens_pop will be continually selected in the next ieration.
## Now select an optimum bandwidth for the variable:  med_niveau_vis 
## The newly selected bandwidth for variable  med_niveau_vis is:  56 
## The bandwidth used in the last ieration is:  56 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  med_niveau_vis seems to be converged and will be kept the same in the following ierations.
## Now select an optimum bandwidth for the variable:  part_log_suroccup 
## The newly selected bandwidth for variable  part_log_suroccup is:  56 
## The bandwidth used in the last ieration is:  56 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  part_log_suroccup seems to be converged and will be kept the same in the following ierations.
## Now select an optimum bandwidth for the variable:  part_agri_nb_emploi 
## The newly selected bandwidth for variable  part_agri_nb_emploi is:  56 
## The bandwidth used in the last ieration is:  56 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  part_agri_nb_emploi seems to be converged and will be kept the same in the following ierations.
## Now select an optimum bandwidth for the variable:  part_cadre_profintellec_nbemploi 
## The newly selected bandwidth for variable  part_cadre_profintellec_nbemploi is:  56 
## The bandwidth used in the last ieration is:  56 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  part_cadre_profintellec_nbemploi seems to be converged and will be kept the same in the following ierations.
##     Ieration  6 the differential change value of RSS (dCVR) is:  0.101719 
## ----------End of    Iteration  6 ----------
##     Iteration  7 :
## Now select an optimum bandwidth for the variable:  Intercept 
## The newly selected bandwidth for variable  Intercept is:  11 
## The bandwidth used in the last ieration is:  11 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  Intercept seems to be converged and will be kept the same in the following ierations.
## Now select an optimum bandwidth for the variable:  perc_maison 
## The newly selected bandwidth for variable  perc_maison is:  56 
## The bandwidth used in the last ieration is:  56 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  perc_maison seems to be converged and will be kept the same in the following ierations.
## Now select an optimum bandwidth for the variable:  dens_pop 
## The newly selected bandwidth for variable  dens_pop is:  41 
## The bandwidth used in the last ieration is:  41 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  dens_pop seems to be converged for  1  times.It will be continually optimized in the next  4  times
##     Ieration  7 the differential change value of RSS (dCVR) is:  0.07537341 
## ----------End of    Iteration  7 ----------
##     Iteration  8 :
## Now select an optimum bandwidth for the variable:  dens_pop 
## The newly selected bandwidth for variable  dens_pop is:  41 
## The bandwidth used in the last ieration is:  41 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  dens_pop seems to be converged for  2  times.It will be continually optimized in the next  3  times
##     Ieration  8 the differential change value of RSS (dCVR) is:  0.06691285 
## ----------End of    Iteration  8 ----------
##     Iteration  9 :
## Now select an optimum bandwidth for the variable:  dens_pop 
## The newly selected bandwidth for variable  dens_pop is:  41 
## The bandwidth used in the last ieration is:  41 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  dens_pop seems to be converged for  3  times.It will be continually optimized in the next  2  times
##     Ieration  9 the differential change value of RSS (dCVR) is:  0.05956372 
## ----------End of    Iteration  9 ----------
##     Iteration  10 :
## Now select an optimum bandwidth for the variable:  dens_pop 
## The newly selected bandwidth for variable  dens_pop is:  41 
## The bandwidth used in the last ieration is:  41 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  dens_pop seems to be converged for  4  times.It will be continually optimized in the next  1  times
##     Ieration  10 the differential change value of RSS (dCVR) is:  0.05253431 
## ----------End of    Iteration  10 ----------
##     Iteration  11 :
## Now select an optimum bandwidth for the variable:  dens_pop 
## The newly selected bandwidth for variable  dens_pop is:  41 
## The bandwidth used in the last ieration is:  41 and the difference between these two bandwidths is:  0 
## The bandwidth for variable  dens_pop seems to be converged and will be kept the same in the following ierations.
##     Ieration  11 the differential change value of RSS (dCVR) is:  0.04631346 
## ----------End of    Iteration  11 ----------
##     Iteration  12 :
##     Ieration  12 the differential change value of RSS (dCVR) is:  0.04089679 
## ----------End of    Iteration  12 ----------
##     Iteration  13 :
##     Ieration  13 the differential change value of RSS (dCVR) is:  0.0361525 
## ----------End of    Iteration  13 ----------
##     Iteration  14 :
##     Ieration  14 the differential change value of RSS (dCVR) is:  0.03197006 
## ----------End of    Iteration  14 ----------
##     Iteration  15 :
##     Ieration  15 the differential change value of RSS (dCVR) is:  0.02826441 
## ----------End of    Iteration  15 ----------
##     Iteration  16 :
##     Ieration  16 the differential change value of RSS (dCVR) is:  0.02496654 
## ----------End of    Iteration  16 ----------
##     Iteration  17 :
##     Ieration  17 the differential change value of RSS (dCVR) is:  0.02201899 
## ----------End of    Iteration  17 ----------
##     Iteration  18 :
##     Ieration  18 the differential change value of RSS (dCVR) is:  0.01937346 
## ----------End of    Iteration  18 ----------
##     Iteration  19 :
##     Ieration  19 the differential change value of RSS (dCVR) is:  0.01698891 
## ----------End of    Iteration  19 ----------
##     Iteration  20 :
##     Ieration  20 the differential change value of RSS (dCVR) is:  0.01483008 
## ----------End of    Iteration  20 ----------
##     Iteration  21 :
##     Ieration  21 the differential change value of RSS (dCVR) is:  0.01286621 
## ----------End of    Iteration  21 ----------
##     Iteration  22 :
##     Ieration  22 the differential change value of RSS (dCVR) is:  0.01106978 
## ----------End of    Iteration  22 ----------
##     Iteration  23 :
##     Ieration  23 the differential change value of RSS (dCVR) is:  0.009415249 
## ----------End of    Iteration  23 ----------
##     Iteration  24 :
##     Ieration  24 the differential change value of RSS (dCVR) is:  0.007877229 
## ----------End of    Iteration  24 ----------
##     Iteration  25 :
##     Ieration  25 the differential change value of RSS (dCVR) is:  0.006427512 
## ----------End of    Iteration  25 ----------
##     Iteration  26 :
##     Ieration  26 the differential change value of RSS (dCVR) is:  0.005028265 
## ----------End of    Iteration  26 ----------
##     Iteration  27 :
##     Ieration  27 the differential change value of RSS (dCVR) is:  0.003610907 
## ----------End of    Iteration  27 ----------
##     Iteration  28 :
##     Ieration  28 the differential change value of RSS (dCVR) is:  0.001957864 
## ----------End of    Iteration  28 ----------
##     Iteration  29 :
##     Ieration  29 the differential change value of RSS (dCVR) is:  0.001713618 
## ----------End of    Iteration  29 ----------
##     Iteration  30 :
##     Ieration  30 the differential change value of RSS (dCVR) is:  0.00278589 
## ----------End of    Iteration  30 ----------
##     Iteration  31 :
##     Ieration  31 the differential change value of RSS (dCVR) is:  0.003322762 
## ----------End of    Iteration  31 ----------
##     Iteration  32 :
##     Ieration  32 the differential change value of RSS (dCVR) is:  0.003619623 
## ----------End of    Iteration  32 ----------
##     Iteration  33 :
##     Ieration  33 the differential change value of RSS (dCVR) is:  0.003769717 
## ----------End of    Iteration  33 ----------
##     Iteration  34 :
##     Ieration  34 the differential change value of RSS (dCVR) is:  0.003819166 
## ----------End of    Iteration  34 ----------
##     Iteration  35 :
##     Ieration  35 the differential change value of RSS (dCVR) is:  0.00379544 
## ----------End of    Iteration  35 ----------
##     Iteration  36 :
##     Ieration  36 the differential change value of RSS (dCVR) is:  0.003716656 
## ----------End of    Iteration  36 ----------
##     Iteration  37 :
##     Ieration  37 the differential change value of RSS (dCVR) is:  0.003595539 
## ----------End of    Iteration  37 ----------
##     Iteration  38 :
##     Ieration  38 the differential change value of RSS (dCVR) is:  0.00344137 
## ----------End of    Iteration  38 ----------
##     Iteration  39 :
##     Ieration  39 the differential change value of RSS (dCVR) is:  0.003261068 
## ----------End of    Iteration  39 ----------
##     Iteration  40 :
##     Ieration  40 the differential change value of RSS (dCVR) is:  0.003059802 
## ----------End of    Iteration  40 ----------
##     Iteration  41 :
##     Ieration  41 the differential change value of RSS (dCVR) is:  0.002841355 
## ----------End of    Iteration  41 ----------
##     Iteration  42 :
##     Ieration  42 the differential change value of RSS (dCVR) is:  0.002608304 
## ----------End of    Iteration  42 ----------
##     Iteration  43 :
##     Ieration  43 the differential change value of RSS (dCVR) is:  0.002362035 
## ----------End of    Iteration  43 ----------
##     Iteration  44 :
##     Ieration  44 the differential change value of RSS (dCVR) is:  0.002102565 
## ----------End of    Iteration  44 ----------
##     Iteration  45 :
##     Ieration  45 the differential change value of RSS (dCVR) is:  0.001827987 
## ----------End of    Iteration  45 ----------
##     Iteration  46 :
##     Ieration  46 the differential change value of RSS (dCVR) is:  0.001533024 
## ----------End of    Iteration  46 ----------
##     Iteration  47 :
##     Ieration  47 the differential change value of RSS (dCVR) is:  0.001204747 
## ----------End of    Iteration  47 ----------
##     Iteration  48 :
##     Ieration  48 the differential change value of RSS (dCVR) is:  0.0008043307 
## ----------End of    Iteration  48 ----------
##     Iteration  49 :
##     Ieration  49 the differential change value of RSS (dCVR) is:  0.000253014 
## ----------End of    Iteration  49 ----------
##     Iteration  50 :
##     Ieration  50 the differential change value of RSS (dCVR) is:  0.0008271053 
## ----------End of    Iteration  50 ----------
##     Iteration  51 :
##     Ieration  51 the differential change value of RSS (dCVR) is:  0.001103321 
## ----------End of    Iteration  51 ----------
##     Iteration  52 :
##     Ieration  52 the differential change value of RSS (dCVR) is:  0.001291758 
## ----------End of    Iteration  52 ----------
##     Iteration  53 :
##     Ieration  53 the differential change value of RSS (dCVR) is:  0.001429585 
## ----------End of    Iteration  53 ----------
##     Iteration  54 :
##     Ieration  54 the differential change value of RSS (dCVR) is:  0.001532507 
## ----------End of    Iteration  54 ----------
##     Iteration  55 :
##     Ieration  55 the differential change value of RSS (dCVR) is:  0.001609177 
## ----------End of    Iteration  55 ----------
##     Iteration  56 :
##     Ieration  56 the differential change value of RSS (dCVR) is:  0.001665106 
## ----------End of    Iteration  56 ----------
##     Iteration  57 :
##     Ieration  57 the differential change value of RSS (dCVR) is:  0.001704137 
## ----------End of    Iteration  57 ----------
##     Iteration  58 :
##     Ieration  58 the differential change value of RSS (dCVR) is:  0.001729128 
## ----------End of    Iteration  58 ----------
##     Iteration  59 :
##     Ieration  59 the differential change value of RSS (dCVR) is:  0.001742305 
## ----------End of    Iteration  59 ----------
##     Iteration  60 :
##     Ieration  60 the differential change value of RSS (dCVR) is:  0.001745454 
## ----------End of    Iteration  60 ----------
##     Iteration  61 :
##     Ieration  61 the differential change value of RSS (dCVR) is:  0.001740053 
## ----------End of    Iteration  61 ----------
##     Iteration  62 :
##     Ieration  62 the differential change value of RSS (dCVR) is:  0.001727341 
## ----------End of    Iteration  62 ----------
##     Iteration  63 :
##     Ieration  63 the differential change value of RSS (dCVR) is:  0.001708376 
## ----------End of    Iteration  63 ----------
##     Iteration  64 :
##     Ieration  64 the differential change value of RSS (dCVR) is:  0.00168407 
## ----------End of    Iteration  64 ----------
##     Iteration  65 :
##     Ieration  65 the differential change value of RSS (dCVR) is:  0.001655218 
## ----------End of    Iteration  65 ----------
##     Iteration  66 :
##     Ieration  66 the differential change value of RSS (dCVR) is:  0.001622514 
## ----------End of    Iteration  66 ----------
##     Iteration  67 :
##     Ieration  67 the differential change value of RSS (dCVR) is:  0.001586571 
## ----------End of    Iteration  67 ----------
##     Iteration  68 :
##     Ieration  68 the differential change value of RSS (dCVR) is:  0.001547928 
## ----------End of    Iteration  68 ----------
##     Iteration  69 :
##     Ieration  69 the differential change value of RSS (dCVR) is:  0.001507062 
## ----------End of    Iteration  69 ----------
##     Iteration  70 :
##     Ieration  70 the differential change value of RSS (dCVR) is:  0.0014644 
## ----------End of    Iteration  70 ----------
##     Iteration  71 :
##     Ieration  71 the differential change value of RSS (dCVR) is:  0.001420316 
## ----------End of    Iteration  71 ----------
##     Iteration  72 :
##     Ieration  72 the differential change value of RSS (dCVR) is:  0.001375146 
## ----------End of    Iteration  72 ----------
##     Iteration  73 :
##     Ieration  73 the differential change value of RSS (dCVR) is:  0.001329187 
## ----------End of    Iteration  73 ----------
##     Iteration  74 :
##     Ieration  74 the differential change value of RSS (dCVR) is:  0.001282702 
## ----------End of    Iteration  74 ----------
##     Iteration  75 :
##     Ieration  75 the differential change value of RSS (dCVR) is:  0.001235925 
## ----------End of    Iteration  75 ----------
##     Iteration  76 :
##     Ieration  76 the differential change value of RSS (dCVR) is:  0.001189062 
## ----------End of    Iteration  76 ----------
##     Iteration  77 :
##     Ieration  77 the differential change value of RSS (dCVR) is:  0.001142296 
## ----------End of    Iteration  77 ----------
##     Iteration  78 :
##     Ieration  78 the differential change value of RSS (dCVR) is:  0.001095787 
## ----------End of    Iteration  78 ----------
##     Iteration  79 :
##     Ieration  79 the differential change value of RSS (dCVR) is:  0.001049675 
## ----------End of    Iteration  79 ----------
##     Iteration  80 :
##     Ieration  80 the differential change value of RSS (dCVR) is:  0.001004084 
## ----------End of    Iteration  80 ----------
##     Iteration  81 :
##     Ieration  81 the differential change value of RSS (dCVR) is:  0.0009591194 
## ----------End of    Iteration  81 ----------
##     Iteration  82 :
##     Ieration  82 the differential change value of RSS (dCVR) is:  0.000914873 
## ----------End of    Iteration  82 ----------
##     Iteration  83 :
##     Ieration  83 the differential change value of RSS (dCVR) is:  0.0008714228 
## ----------End of    Iteration  83 ----------
##     Iteration  84 :
##     Ieration  84 the differential change value of RSS (dCVR) is:  0.0008288348 
## ----------End of    Iteration  84 ----------
##     Iteration  85 :
##     Ieration  85 the differential change value of RSS (dCVR) is:  0.0007871638 
## ----------End of    Iteration  85 ----------
##     Iteration  86 :
##     Ieration  86 the differential change value of RSS (dCVR) is:  0.0007464546 
## ----------End of    Iteration  86 ----------
##     Iteration  87 :
##     Ieration  87 the differential change value of RSS (dCVR) is:  0.0007067425 
## ----------End of    Iteration  87 ----------
##     Iteration  88 :
##     Ieration  88 the differential change value of RSS (dCVR) is:  0.0006680546 
## ----------End of    Iteration  88 ----------
##     Iteration  89 :
##     Ieration  89 the differential change value of RSS (dCVR) is:  0.0006304099 
## ----------End of    Iteration  89 ----------
##     Iteration  90 :
##     Ieration  90 the differential change value of RSS (dCVR) is:  0.0005938205 
## ----------End of    Iteration  90 ----------
##     Iteration  91 :
##     Ieration  91 the differential change value of RSS (dCVR) is:  0.0005582915 
## ----------End of    Iteration  91 ----------
##     Iteration  92 :
##     Ieration  92 the differential change value of RSS (dCVR) is:  0.0005238218 
## ----------End of    Iteration  92 ----------
##     Iteration  93 :
##     Ieration  93 the differential change value of RSS (dCVR) is:  0.0004904043 
## ----------End of    Iteration  93 ----------
##     Iteration  94 :
##     Ieration  94 the differential change value of RSS (dCVR) is:  0.0004580259 
## ----------End of    Iteration  94 ----------
##     Iteration  95 :
##     Ieration  95 the differential change value of RSS (dCVR) is:  0.0004266676 
## ----------End of    Iteration  95 ----------
##     Iteration  96 :
##     Ieration  96 the differential change value of RSS (dCVR) is:  0.0003963045 
## ----------End of    Iteration  96 ----------
##     Iteration  97 :
##     Ieration  97 the differential change value of RSS (dCVR) is:  0.0003669051 
## ----------End of    Iteration  97 ----------
##     Iteration  98 :
##     Ieration  98 the differential change value of RSS (dCVR) is:  0.0003384308 
## ----------End of    Iteration  98 ----------
##     Iteration  99 :
##     Ieration  99 the differential change value of RSS (dCVR) is:  0.0003108344 
## ----------End of    Iteration  99 ----------
##     Iteration  100 :
##     Ieration  100 the differential change value of RSS (dCVR) is:  0.0002840582 
## ----------End of    Iteration  100 ----------
##     Iteration  101 :
##     Ieration  101 the differential change value of RSS (dCVR) is:  0.0002580309 
## ----------End of    Iteration  101 ----------
##     Iteration  102 :
##     Ieration  102 the differential change value of RSS (dCVR) is:  0.0002326619 
## ----------End of    Iteration  102 ----------
##     Iteration  103 :
##     Ieration  103 the differential change value of RSS (dCVR) is:  0.0002078331 
## ----------End of    Iteration  103 ----------
##     Iteration  104 :
##     Ieration  104 the differential change value of RSS (dCVR) is:  0.0001833824 
## ----------End of    Iteration  104 ----------
##     Iteration  105 :
##     Ieration  105 the differential change value of RSS (dCVR) is:  0.0001590743 
## ----------End of    Iteration  105 ----------
##     Iteration  106 :
##     Ieration  106 the differential change value of RSS (dCVR) is:  0.0001345353 
## ----------End of    Iteration  106 ----------
##     Iteration  107 :
##     Ieration  107 the differential change value of RSS (dCVR) is:  0.0001090957 
## ----------End of    Iteration  107 ----------
##     Iteration  108 :
##     Ieration  108 the differential change value of RSS (dCVR) is:  8.128309e-05 
## ----------End of    Iteration  108 ----------
##     Iteration  109 :
##     Ieration  109 the differential change value of RSS (dCVR) is:  4.60327e-05 
## ----------End of    Iteration  109 ----------
##     Iteration  110 :
##     Ieration  110 the differential change value of RSS (dCVR) is:  4.062144e-05 
## ----------End of    Iteration  110 ----------
##     Iteration  111 :
##     Ieration  111 the differential change value of RSS (dCVR) is:  6.914646e-05 
## ----------End of    Iteration  111 ----------
##     Iteration  112 :
##     Ieration  112 the differential change value of RSS (dCVR) is:  8.5722e-05 
## ----------End of    Iteration  112 ----------
##     Iteration  113 :
##     Ieration  113 the differential change value of RSS (dCVR) is:  9.70491e-05 
## ----------End of    Iteration  113 ----------
##     Iteration  114 :
##     Ieration  114 the differential change value of RSS (dCVR) is:  0.0001051343 
## ----------End of    Iteration  114 ----------
##     Iteration  115 :
##     Ieration  115 the differential change value of RSS (dCVR) is:  0.0001109407 
## ----------End of    Iteration  115 ----------
##     Iteration  116 :
##     Ieration  116 the differential change value of RSS (dCVR) is:  0.0001150325 
## ----------End of    Iteration  116 ----------
##     Iteration  117 :
##     Ieration  117 the differential change value of RSS (dCVR) is:  0.0001177807 
## ----------End of    Iteration  117 ----------
##     Iteration  118 :
##     Ieration  118 the differential change value of RSS (dCVR) is:  0.000119448 
## ----------End of    Iteration  118 ----------
##     Iteration  119 :
##     Ieration  119 the differential change value of RSS (dCVR) is:  0.0001202308 
## ----------End of    Iteration  119 ----------
##     Iteration  120 :
##     Ieration  120 the differential change value of RSS (dCVR) is:  0.000120281 
## ----------End of    Iteration  120 ----------
##     Iteration  121 :
##     Ieration  121 the differential change value of RSS (dCVR) is:  0.0001197199 
## ----------End of    Iteration  121 ----------
##     Iteration  122 :
##     Ieration  122 the differential change value of RSS (dCVR) is:  0.0001186466 
## ----------End of    Iteration  122 ----------
##     Iteration  123 :
##     Ieration  123 the differential change value of RSS (dCVR) is:  0.0001171432 
## ----------End of    Iteration  123 ----------
##     Iteration  124 :
##     Ieration  124 the differential change value of RSS (dCVR) is:  0.0001152787 
## ----------End of    Iteration  124 ----------
##     Iteration  125 :
##     Ieration  125 the differential change value of RSS (dCVR) is:  0.0001131121 
## ----------End of    Iteration  125 ----------
##     Iteration  126 :
##     Ieration  126 the differential change value of RSS (dCVR) is:  0.0001106937 
## ----------End of    Iteration  126 ----------
##     Iteration  127 :
##     Ieration  127 the differential change value of RSS (dCVR) is:  0.000108067 
## ----------End of    Iteration  127 ----------
##     Iteration  128 :
##     Ieration  128 the differential change value of RSS (dCVR) is:  0.0001052699 
## ----------End of    Iteration  128 ----------
##     Iteration  129 :
##     Ieration  129 the differential change value of RSS (dCVR) is:  0.0001023354 
## ----------End of    Iteration  129 ----------
##     Iteration  130 :
##     Ieration  130 the differential change value of RSS (dCVR) is:  9.929239e-05 
## ----------End of    Iteration  130 ----------
##     Iteration  131 :
##     Ieration  131 the differential change value of RSS (dCVR) is:  9.616616e-05 
## ----------End of    Iteration  131 ----------
##     Iteration  132 :
##     Ieration  132 the differential change value of RSS (dCVR) is:  9.297901e-05 
## ----------End of    Iteration  132 ----------
##     Iteration  133 :
##     Ieration  133 the differential change value of RSS (dCVR) is:  8.975046e-05 
## ----------End of    Iteration  133 ----------
##     Iteration  134 :
##     Ieration  134 the differential change value of RSS (dCVR) is:  8.649773e-05 
## ----------End of    Iteration  134 ----------
##     Iteration  135 :
##     Ieration  135 the differential change value of RSS (dCVR) is:  8.323595e-05 
## ----------End of    Iteration  135 ----------
##     Iteration  136 :
##     Ieration  136 the differential change value of RSS (dCVR) is:  7.997836e-05 
## ----------End of    Iteration  136 ----------
##     Iteration  137 :
##     Ieration  137 the differential change value of RSS (dCVR) is:  7.673659e-05 
## ----------End of    Iteration  137 ----------
##     Iteration  138 :
##     Ieration  138 the differential change value of RSS (dCVR) is:  7.352076e-05 
## ----------End of    Iteration  138 ----------
##     Iteration  139 :
##     Ieration  139 the differential change value of RSS (dCVR) is:  7.033971e-05 
## ----------End of    Iteration  139 ----------
##     Iteration  140 :
##     Ieration  140 the differential change value of RSS (dCVR) is:  6.720106e-05 
## ----------End of    Iteration  140 ----------
##     Iteration  141 :
##     Ieration  141 the differential change value of RSS (dCVR) is:  6.411143e-05 
## ----------End of    Iteration  141 ----------
##     Iteration  142 :
##     Ieration  142 the differential change value of RSS (dCVR) is:  6.107632e-05 
## ----------End of    Iteration  142 ----------
##     Iteration  143 :
##     Ieration  143 the differential change value of RSS (dCVR) is:  5.810054e-05 
## ----------End of    Iteration  143 ----------
##     Iteration  144 :
##     Ieration  144 the differential change value of RSS (dCVR) is:  5.518801e-05 
## ----------End of    Iteration  144 ----------
##     Iteration  145 :
##     Ieration  145 the differential change value of RSS (dCVR) is:  5.234191e-05 
## ----------End of    Iteration  145 ----------
##     Iteration  146 :
##     Ieration  146 the differential change value of RSS (dCVR) is:  4.956483e-05 
## ----------End of    Iteration  146 ----------
##     Iteration  147 :
##     Ieration  147 the differential change value of RSS (dCVR) is:  4.685875e-05 
## ----------End of    Iteration  147 ----------
##     Iteration  148 :
##     Ieration  148 the differential change value of RSS (dCVR) is:  4.422502e-05 
## ----------End of    Iteration  148 ----------
##     Iteration  149 :
##     Ieration  149 the differential change value of RSS (dCVR) is:  4.166458e-05 
## ----------End of    Iteration  149 ----------
##     Iteration  150 :
##     Ieration  150 the differential change value of RSS (dCVR) is:  3.917786e-05 
## ----------End of    Iteration  150 ----------
##     Iteration  151 :
##     Ieration  151 the differential change value of RSS (dCVR) is:  3.676487e-05 
## ----------End of    Iteration  151 ----------
##     Iteration  152 :
##     Ieration  152 the differential change value of RSS (dCVR) is:  3.442506e-05 
## ----------End of    Iteration  152 ----------
##     Iteration  153 :
##     Ieration  153 the differential change value of RSS (dCVR) is:  3.215768e-05 
## ----------End of    Iteration  153 ----------
##     Iteration  154 :
##     Ieration  154 the differential change value of RSS (dCVR) is:  2.996143e-05 
## ----------End of    Iteration  154 ----------
##     Iteration  155 :
##     Ieration  155 the differential change value of RSS (dCVR) is:  2.783453e-05 
## ----------End of    Iteration  155 ----------
##     Iteration  156 :
##     Ieration  156 the differential change value of RSS (dCVR) is:  2.577497e-05 
## ----------End of    Iteration  156 ----------
##     Iteration  157 :
##     Ieration  157 the differential change value of RSS (dCVR) is:  2.377999e-05 
## ----------End of    Iteration  157 ----------
##     Iteration  158 :
##     Ieration  158 the differential change value of RSS (dCVR) is:  2.184639e-05 
## ----------End of    Iteration  158 ----------
##     Iteration  159 :
##     Ieration  159 the differential change value of RSS (dCVR) is:  1.997021e-05 
## ----------End of    Iteration  159 ----------
##     Iteration  160 :
##     Ieration  160 the differential change value of RSS (dCVR) is:  1.814646e-05 
## ----------End of    Iteration  160 ----------
##     Iteration  161 :
##     Ieration  161 the differential change value of RSS (dCVR) is:  1.636901e-05 
## ----------End of    Iteration  161 ----------
##     Iteration  162 :
##     Ieration  162 the differential change value of RSS (dCVR) is:  1.462952e-05 
## ----------End of    Iteration  162 ----------
##     Iteration  163 :
##     Ieration  163 the differential change value of RSS (dCVR) is:  1.291706e-05 
## ----------End of    Iteration  163 ----------
##     Iteration  164 :
##     Ieration  164 the differential change value of RSS (dCVR) is:  1.121533e-05 
## ----------End of    Iteration  164 ----------
##     Iteration  165 :
##     Ieration  165 the differential change value of RSS (dCVR) is:  9.498759e-06 
## ----------End of    Iteration  165 ----------

mgwr.bw  <- round(MGWR[[2]]$bws,1) # Nombre de voisins pour chaque prédicteur
#mgwr.bw

# Exploration des résultats statistiques
print(MGWR)

##    ***********************************************************************
##    *                       Package   GWmodel                             *
##    ***********************************************************************
##    Program starts at: 2023-10-07 14:27:54.994629 
##    Call:
##    gwr.multiscale(formula = prix_med ~ perc_log_vac + perc_maison + 
##     perc_tiny_log + dens_pop + med_niveau_vis + part_log_suroccup + 
##     part_agri_nb_emploi + part_cadre_profintellec_nbemploi, data = data_immo_bzh_sp, 
##     kernel = "gaussian", adaptive = TRUE, bws0 = rep(1, 4), predictor.centered = rep(T, 
##         3))
## 
##    Dependent (y) variable:  prix_med
##    Independent variables:  perc_log_vac perc_maison perc_tiny_log dens_pop med_niveau_vis part_log_suroccup part_agri_nb_emploi part_cadre_profintellec_nbemploi
##    Number of data points: 57
##    ***********************************************************************
##    *                       Multiscale (PSDM) GWR                          *
##    ***********************************************************************
## 
##    *********************Model calibration information*********************
##    Kernel function: gaussian 
##    Adaptive bandwidths for each coefficient(number of nearest neighbours): 
##               (Intercept) perc_log_vac perc_maison perc_tiny_log dens_pop
##    Bandwidth           11           53          56            39       41
##               med_niveau_vis part_log_suroccup part_agri_nb_emploi
##    Bandwidth              56                56                  56
##               part_cadre_profintellec_nbemploi
##    Bandwidth                                56
## 
##    ****************Summary of GWR coefficient estimates:******************
##                                         Min.  1st Qu.   Median  3rd Qu.
##    Intercept                        1649.395 1723.184 1865.122 1929.827
##    perc_log_vac                     -302.118 -298.338 -296.688 -294.292
##    perc_maison                      -591.452 -591.063 -588.224 -583.696
##    perc_tiny_log                    -413.370 -408.776 -402.613 -395.289
##    dens_pop                         -219.688 -217.161 -191.024 -167.340
##    med_niveau_vis                    339.684  344.076  345.308  347.536
##    part_log_suroccup                 339.851  340.086  342.984  346.007
##    part_agri_nb_emploi                98.592  102.212  103.293  104.229
##    part_cadre_profintellec_nbemploi  -64.223  -63.736  -59.797  -55.391
##                                         Max.
##    Intercept                        1952.455
##    perc_log_vac                     -292.470
##    perc_maison                      -582.969
##    perc_tiny_log                    -392.117
##    dens_pop                         -163.384
##    med_niveau_vis                    353.134
##    part_log_suroccup                 346.272
##    part_agri_nb_emploi               106.722
##    part_cadre_profintellec_nbemploi  -54.977
##    ************************Diagnostic information*************************
##    Number of data points: 57 
##    Effective number of parameters (2trace(S) - trace(S'S)): 14.16268 
##    Effective degrees of freedom (n-2trace(S) + trace(S'S)): 42.83732 
##    AICc value:  788.7976 
##    AIC value:  765.7213 
##    BIC value:  745.9901 
##    Residual sum of squares:  1837688 
##    R-square value:  0.8802325 
##    Adjusted R-square value:  0.839689 
## 
##    ***********************************************************************
##    Program stops at: 2023-10-07 14:27:58.977106

Il me semble que le nombre de plus proches voisins (number of nearest neighbours) nous indique si l’effet de notre prédicteurs relève d’un processus qui soit plus global (càd échelle du territoire) ou au contraire beaucoup plus localisé. Plus le nombre de voisins est petit plus l’effet est localisé.

Pour visualiser les Betas :

datatable(as.data.frame(MGWR$SDF))

A partir de ces résultats, on peut refaire toutes les analyses et cartes réalisées avec la GWR standard.

Par exemple, les cartes des coefficients des VI par epci Breton

data_immo_bzh$agri.mgwr=MGWR$SDF$part_agri_nb_emploi
data_immo_bzh$perc_maison.mgwr <- MGWR$SDF$perc_maison
data_immo_bzh$dens_pop.mgwr=MGWR$SDF$dens_pop
data_immo_bzh$med_vie.mgwr=MGWR$SDF$med_niveau_vis
data_immo_bzh$logvac.mgwr=MGWR$SDF$perc_log_vac
data_immo_bzh$tinylog.mgwr=MGWR$SDF$perc_tiny_log
data_immo_bzh$suroccup.mgwr=MGWR$SDF$part_log_suroccup
data_immo_bzh$cadre.mgwr=MGWR$SDF$part_cadre_profintellec_nbemploi

# Réaliser la collection des cartes

par(mfrow = c(2, 4))

mf_map(x = data_immo_bzh, var = "agri.mgwr", type = "choro", pal= "Earth")
mf_title("Coefficients Agriculteurs MGWR")

mf_map(x = data_immo_bzh, var = "perc_maison.mgwr", type = "choro", pal= "Earth")
mf_title("Coefficients de Maison MGWR")

mf_map(x = data_immo_bzh, var = "dens_pop.mgwr", type = "choro", pal= "Earth")
mf_title("Coefficients de dens pop MGWR")

mf_map(x = data_immo_bzh, var = "med_vie.mgwr", type = "choro", pal= "Earth")
mf_title("Coefficients de Médianne niveau de vie MGWR")

mf_map(x = data_immo_bzh, var = "logvac.mgwr", type = "choro", pal= "Earth")
mf_title("Coefficients de Logements vacants MGWR")

mf_map(x = data_immo_bzh, var = "tinylog.mgwr", type = "choro", pal= "Earth")
mf_title("Coefficients de Petits logements MGWR")

mf_map(x = data_immo_bzh, var = "suroccup.mgwr", type = "choro", pal= "Earth")
mf_title("Coefficients de logement suroocupés MGWR")

mf_map(x = data_immo_bzh, var = "cadre.mgwr", type = "choro", pal= "Earth")
mf_title("Coefficients de Cadre MGWR")

3.2 Régionalisation

Ici l’objectif va être de définir des sous-espaces sur la base des coefficients de la GWR, sous-espaces qui vont se caractériser par l’homogénéité des coefficients des prédicteurs au sein du sous-espace.

Thierry feuillet décrit cette méthode de cette manière :

“Ce processus de découpage de l’espace en sous-régions homogènes se nomme la régionalisation. C’est une extension de la classification classique : on y ajoute un critère de contiguité spatiale. La régionalisation est donc une classification spatiale.

Il existe plusieurs méthodes de régionalisation. Un des principes les plus répandus consiste à établir une classification à la fois sur la base de la ressemblance entre les observations, et sur leur proximité dans l’espace géographique.

Nous allons ici utiliser l’algorithme SKATER (Spatial Klustering Analysis by Tree Edge Removal), méthode proposée par Assunçao et al. (2006) et déjà appliqué dans un contexte de recherche similaire au notre par Helbich et al. (2013). Par ailleurs, une description très pédagogique de la méthode est disponible ici : http://www.jms-insee.fr/2018/S08_5_ACTE_ROUSSEZ_JMS2018.pdf

L’algorithme SKATER comporte 4 étapes (cf. doc cité ci-dessus) :

1- Constuction d’un graphe de voisinage (contiguité ou knn) 2- Pondération des liens du graphe à partir de la matrice de dissimilarité 3- Construction de l’arbre portant minimal, en retenant le lien avec le voisin le plus ressemblant pour chaque noeud 4- Elagage de l’arbre maximisant la variance inter-classes des sous-graphes

Avant de poursuivre on va donc recalculer un modèle GWR pour notre exemple breton:

library(spgwr)

## NOTE: This package does not constitute approval of GWR
## as a method of spatial analysis; see example(gwr)

best.bzh <- gwr.sel(formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + part_cadre_profintellec_nbemploi, 
                    data = data_immo_bzh, 
                    coords = st_coordinates(st_centroid(data_immo_bzh)))

## Warning: st_centroid assumes attributes are constant over geometries

## Bandwidth: 109909.9 CV score: 4736532 
## Bandwidth: 177660.5 CV score: 4681812 
## Bandwidth: 219532.6 CV score: 4711813 
## Bandwidth: 172830.3 CV score: 4678660 
## Bandwidth: 148796.9 CV score: 4668839 
## Bandwidth: 136640.9 CV score: 4671999 
## Bandwidth: 149754.5 CV score: 4668918 
## Bandwidth: 147704.4 CV score: 4668797 
## Bandwidth: 147354.2 CV score: 4668795 
## Bandwidth: 147390.6 CV score: 4668795 
## Bandwidth: 147393.9 CV score: 4668795 
## Bandwidth: 147393.8 CV score: 4668795 
## Bandwidth: 147393.8 CV score: 4668795 
## Bandwidth: 147392.5 CV score: 4668795 
## Bandwidth: 147393.3 CV score: 4668795 
## Bandwidth: 147393.6 CV score: 4668795 
## Bandwidth: 147393.7 CV score: 4668795 
## Bandwidth: 147393.7 CV score: 4668795 
## Bandwidth: 147393.7 CV score: 4668795 
## Bandwidth: 147393.7 CV score: 4668795 
## Bandwidth: 147393.7 CV score: 4668795 
## Bandwidth: 147393.7 CV score: 4668795

bzh_gwr <- gwr(formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + part_cadre_profintellec_nbemploi, 
               data = data_immo_bzh, 
               coord = st_coordinates(st_centroid(data_immo_bzh)),
               bandwidth = best.bzh, 
               gweight = gwr.Gauss,
               hatmatrix = TRUE)

## Warning: st_centroid assumes attributes are constant over geometries

3.2.1 Première étape : préparation de la table des coefficients GWR

On fait une standardisation (centrage-réduction) pour rendre nos différents coefficients comparables.

bzh.stand <- bzh_gwr$SDF %>% 
  as.data.frame() %>% 
  select(3:10) %>% # On ne conserve que les colonne des coeff de nos prédicteurs
  mutate_all(~scale(.)) %>% 
  rename_with(~paste(.x, "b", sep = "_"))

data_bzh <- cbind(data_immo_bzh, bzh.stand)

data_bzh_sp <- as(data_bzh, "Spatial")

3.2.2 Computation de l’algorithme SKATER : Définition du voisinage de chaque point

De ce que j’ai compris, Skater est un algorithme pour faire du clustering de données spatiales, en s’assurant que les clusters sont constitués d’objets contigus. Pour en savoir plus sur skater, voir ce post (skater?).

knn <- knearneigh(bzh_gwr$SDF, k = 50) # On utilise knearneigh car nous ne sommes plus sur un shape avec des polygons - on avait alors utilisé la fonction poly2nb - mais sur une matrice de points coordonnées

## Warning in knearneigh(bzh_gwr$SDF, k = 50): k greater than one-third of the
## number of data points

nb <- knn2nb(knn)
plot(nb, coords = coordinates(bzh_gwr$SDF), col="blue")

Calibrage du coût des arêtes et de la pondération spatiale

costs <- nbcosts(nb, data = bzh.stand)
costsW <- nb2listw(nb, costs, style="B")

Minimisation de l’arbre et classification

costsTree <- mstree(costsW)
plot(costsTree, coords = coordinates(bzh_gwr$SDF), col="blue", main = "Arbre portant minimal")

Ici définition en 5 clusters, c’est arbitraire mais orienté par les analyses précédentes.

# Mettre spdep:: devant la fonction car le package rgeoda possède la même fonction qui fait la même chose mais pas les mêmes arguments
clus5 <- spdep::skater(edges = costsTree[,1:2], data = bzh.stand, ncuts = 5)

bzhClus <- data_immo_bzh %>% 
mutate(clus = as.factor(clus5$groups)) %>% 
bind_cols(bzh.stand)

mf_map(x = bzhClus, var = "clus", type = "typo", pal= "Set 2")
mf_title("régionalisation Bretagne")

Et grâce au code qui suit vous pouvez caractériser les clusters.

library(ggplot2)

nomVar <- c("perc_log_vac_b","perc_maison_b","perc_tiny_log_b","dens_pop_b","med_niveau_vis_b", "part_log_suroccup_b","part_agri_nb_emploi_b","part_cadre_profintellec_nbemploi_b")

clusProfile <- bzhClus[, c(nomVar, "clus")] %>% 
  group_by(clus) %>% 
  summarise_each(funs(mean)) %>% 
  st_drop_geometry()

## Warning: `summarise_each_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `across()` instead.
## ℹ The deprecated feature was likely used in the dplyr package.
##   Please report the issue at <https://github.com/tidyverse/dplyr/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `funs()` was deprecated in dplyr 0.8.0.
## ℹ Please use a list of either functions or lambdas:
## 
## # Simple named list: list(mean = mean, median = median)
## 
## # Auto named with `tibble::lst()`: tibble::lst(mean, median)
## 
## # Using lambdas list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: There were 6 warnings in `summarise()`.
## The first warning was:
## ℹ In argument: `geometry = mean(geometry)`.
## ℹ In group 1: `clus = 1`.
## Caused by warning in `mean.default()`:
## ! l'argument n'est ni numérique, ni logique : renvoi de NA
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 5 remaining warnings.

clusLong <- reshape2::melt(clusProfile, id.vars = "clus")

profilePlot <- ggplot(clusLong) +
  geom_bar(aes(x = variable, y = value), 
           fill = "grey25", 
           position = "identity", 
           stat = "identity") +
  scale_x_discrete("Effet") + 
  scale_y_continuous("Valeur moyenne par classe") +
  facet_wrap(~ clus) + 
  coord_flip() + 
  theme(strip.background = element_rect(fill="grey25"),
        strip.text = element_text(colour = "grey85", face = "bold"))

profilePlot

Bibliographie

ANSELIN, Luc, IBNU, Syabri et YOUNGIHN, Kho, 2006. GeoDa: An Introduction to Spatial Data Analysis. [en ligne]. 2006. S.l. : s.n. Disponible à l'adresse : https://geodacenter.github.io/.

BAILEY, Trevor C. et GATRELL, Anthony C., 1995. Interactive spatial data analysis. S.l. : s.n.

BRUNSDON, Chris, FOTHERINGHAM, A. Stewart et CHARLTON, Martin E., 1996. Geographically Weighted Regression: A Method for Exploring Spatial Nonstationarity. In : Journal of the Royal Statistical Society: Series D (The Statistician) [en ligne]. 1996. Vol. 28, n° 4. DOI 10.1111/j.1538-4632.1996.tb00936.x. Disponible à l'adresse : https://doi.org/10.1111/j.1538-4632.1996.tb00936.x.

DE BELLEFON, Marie-Pierre, LOONIS, Vincent et LE GLEUT, Renan, 2018. Codifier la structure de voisinage. In : EUROSTAT, Insee - (éd.), Manuel d’analyse spatiale [en ligne]. S.l. : s.n. pp. 33‑50. Disponible à l'adresse : https://www.insee.fr/fr/information/3635442.

FEUILLET, Thierry, COSSART, Étienne et COMMENGES, Hadrien, 2019. Manuel de géographie quantitative - Concepts, outils, méthodes. S.l. : s.n.

LE CAMPION, Grégoire, 2021. Analyse des corrélations avec easystats. In : RZine [en ligne]. 2021. DOI 10.48645/QHAV-CB52. Disponible à l'adresse : https://rzine.fr/publication_rzine/20201127_glecampion_initiation_aux_correlations/.

MATHIAN, Hélène et PIRON, Marie, 2001. Échelles géographiques et méthodes statistiques multidimensionnelles. In : SANDERS, Lena (éd.), Modèles en analyse spatiale. S.l. : Hermès - Lavoisier. pp. 61‑103.

OLIVEAU, Sébastien, 2011. L’espace compte ! Mesurer les structures spatiales du changement social. Habilitation à diriger des recherches. S.l. : Université Aix-Marseille 1, UFR des sciences géographiques et de l’aménagement.

Annexes

Info session

setting	value
version	R version 4.3.1 (2023-06-16)
os	Ubuntu 22.04.3 LTS
system	x86_64, linux-gnu
ui	X11
language	fr_FR:en
collate	fr_FR.UTF-8
ctype	fr_FR.UTF-8
tz	Europe/Paris
date	2023-10-07
pandoc	3.1.1 @ /usr/lib/rstudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)

package	ondiskversion	source
car	3.1.2	CRAN (R 4.3.0)
carData	3.0.5	CRAN (R 4.3.0)
correlation	0.8.4	CRAN (R 4.3.0)
corrplot	0.92	CRAN (R 4.3.0)
digest	0.6.33	CRAN (R 4.3.1)
dplyr	1.1.2	CRAN (R 4.3.0)
DT	0.28	CRAN (R 4.3.0)
GGally	2.1.2	CRAN (R 4.3.0)
ggplot2	3.4.2	CRAN (R 4.3.0)
gtsummary	1.7.1	CRAN (R 4.3.0)
GWmodel	2.2.9	CRAN (R 4.3.0)
here	1.0.1	CRAN (R 4.3.0)
mapsf	0.6.1	CRAN (R 4.3.0)
maptools	1.1.7	CRAN (R 4.3.0)
Matrix	1.6.0	CRAN (R 4.3.1)
plotly	4.10.2	CRAN (R 4.3.0)
RColorBrewer	1.1.3	CRAN (R 4.3.0)
Rcpp	1.0.11	CRAN (R 4.3.1)
rgeoda	0.0.10.2	CRAN (R 4.3.0)
robustbase	0.99.0	CRAN (R 4.3.0)
sf	1.0.14	CRAN (R 4.3.1)
sp	2.0.0	CRAN (R 4.3.0)
spatialreg	1.2.9	CRAN (R 4.3.0)
spData	2.2.2	CRAN (R 4.3.0)
spdep	1.2.8	CRAN (R 4.3.0)
spgwr	0.6.36	CRAN (R 4.3.1)